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APRESENTAÇÃO 


Este livro foi desenvolvido para dar suporte ao estudo de uma disciplina 
básica de Cálculo Numérico. 

O Cálculo Numérico tem por objetivo o estudo de métodos numéricos para 
a resolução de problemas em geral e, em particular, destacamos neste texto 
os problemas ambientais. Os métodos numéricos serão desenvolvidos visando 
eficiência e estabilidade durante a aplicação em problemas gerais. 

Consideramos de fundamental importância o Cálculo Numérico no que 
se refere à aplicação das disciplinas básicas como Cálculo Diferencial Integral, 
Geometria Analítica, Álgebra Linear, Introdução à Computação, entre outras, e 
a utilização dos conceitos e resultados adquiridos nessas disciplinas, bem como a 
aplicação na resolução numérica de problemas. 

Apresentamos métodos numéricos para solução de sistemas de equações 
lineares, equações algébricas, interpolação e ajuste de funções, integração nu¬ 
mérica e equações diferenciais. 

Detalhamos alguns aspectos teóricos e computacionais, além da utiliza¬ 
ção de softwares especializados e reconhecidos na literatura, com a certeza de 
que estaremos contribuindo com ferramentas importantes para a elaboração 
de projetos e pesquisas que estarão presentes nas atividades dos profissionais de 
ciências em geral. 

Por fim, agradecemos a Carla Taviani Lucke da Silva Ghidini e Silvia Maria 
Pereira Grandi dos Santos pela colaboração na revisão deste livro. 
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UNIDADE 1 


Teoria dos Erros em Processos Numéricos 



1.1 Introdução 


Nas diversas áreas científicas, diante da resolução numérica de um pro¬ 
blema, deparamo-nos com a necessidade de tratar com soluções aproximadas, 
desde a interpretação do problema, a modelagem matemática correspondente, 
erros nos dados obtidos por meio de medidas experimentais e implementação 
de dados no computador. 

Para melhor entendimento, resumimos alguns passos considerados impor¬ 
tantes na resolução de um problema da seguinte forma: 

1. Apresentação do problema real, entendimento e levantamento dos dados; 

2. Formulação matemática correspondente do problema-modelo matemático; 

3. Elaboração de um plano de resolução e a escolha de um método nu¬ 
mérico adequado; 

4. Implementação computacional do método, linguagem de programação; 

5. Análise sobre a coerência dos resultados obtidos e o problema inicial 
proposto; 

6. Reformular o modelo matemático e/ou escolher novo método numérico 
de resolução, caso o Passo 5 não esteja satisfeito. 


A execução dos passos anteriores é chamada de Validação do Modelo 
e as aproximações consideradas nesses passos levam a alguns tipos de erros, 
conforme exemplo a seguir: 


Exemplo 1.1 



Considere o problema de transporte de lixo contaminado, conforme exibido 
na Figura 1: 


Figura 1 Lixo contaminado. 






Uma transportadora possui três tipos de caminhões representados por C,, 
C 2 e C 3 , os quais são equipados para levar três tipos diferentes de materiais con¬ 
taminados M i; M 2 e M 3 para lixões adequados, conforme a Tabela 1.1: 


Tabela 1.1 



Ml 

M2 

M3 

Cl 

6.8 

1.9 

1.7 

C2 

1.5 

2.6 

1.4 

C3 

2.3 

1.8 

2.3 


O caminhão C 1 pode levar 6.8 toneladas do tipo M 15 1.9 tonelada do tipo M 2 
e 1.7 tonelada do M 3 . Deseja-se saber quantos caminhões de cada tipo devemos 
enviar para transportar 16 toneladas do tipo M r 7 toneladas do M 2 e 6 toneladas 
do M 3 ? 

Exibiremos a execução de todos os passos citados anteriormente e fare¬ 
mos uma análise final dos resultados obtidos. 

• Passo 1 

Neste passo, teremos o entendimento do problema observando todos os 
detalhes, características e condições intrínsecas do mesmo. 

• Passo 2 

Faremos a formulação matemática do problema, para isso, devemos deter¬ 
minar as variáveis do problema. 

Assim, seja o número de caminhões x jt j = 1,2,3 a ser determinado para 
o transporte do lixo contaminado. Dessa forma, observando a tabela de dados 
vemos que podemos escrever a seguinte equação para o transporte do lixo con¬ 
taminado 

6.8 x 1 +1.5 x 2 +2.3 x 3 =16 

Seguindo o mesmo raciocínio, para o transporte dos lixos contaminados 
M 2 e M 3 , temos o seguinte sistema de equações lineares para ser resolvido, que 
será estudado nas próximas unidades de estudo. 

6.8x 1 + 1.5 x 2 + 2.3x 3 = 16 
< 1.9x 1 + 2.6 x 2 + 1.8x 3 = 7 
1.7x 1 + 1.4 x 2 + 2.3x 3 = 6 
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Passo 3 


Devemos escolher o método adequado de resolução. Neste caso, vamos 
utilizar o Método de Eliminação de Gauss, o qual será alvo de estudo nas pró¬ 
ximas unidades. 

• Passo 4 

Neste passo, podemos implementar no computador o método de resolução 
escolhido ou utilizar um software apropriado e teremos a solução do sistema 
como segue: 


x 1 = 1.9452, x 2 = 0.7953 e x 3 = 0.6868 


• Passo 5 

Temos que analisar a solução encontrada, pois o número de caminhões 
deverá ser um número inteiro positivo. Como temos uma solução aproximada 
não inteira para o sistema de equações lineares, podemos aproximá-la por um 
número inteiro mais próximo de tal solução. Dessa forma, tomamos como solu¬ 
ção x 1 =2.0, x 2 =1.0 e x 3 = 1.0. 

Assim, devemos contratar dois caminhões do tipo C i; um caminhão do tipo 
C 2 e um caminhão do tipo C 3 . 

Observações 

a) Na modelagem matemática do problema, consideramos uma aproxima¬ 
ção para o problema dado, uma vez que a capacidade de transporte 
de cada material em cada caminhão é estimada, isto é, com uma certa 
margem de erro. 

Alguns detalhes foram desconsiderados, como a integralidade da solução 
e custos envolvidos no processo de transporte. 

b) A solução obtida durante a execução do método de resolução foi apre¬ 
sentada com uma aproximação de quatro casas decimais, o que signifi¬ 
ca que estamos cometendo erros em todas essas aproximações. 

Diante das considerações anteriores sobre o tratamento dos problemas com 
aproximações, desde a modelagem matemática, os erros cometidos durante o 
processamento dos métodos de resolução no computador, é necessário um en¬ 
tendimento geral sobre tipos de erros existentes, como descreveremos a seguir. 




1.2 Erros na fase da modelagem 


Devido às simplificações no processo de modelagem matemática de um 
problema, que muitas vezes são necessárias, podem ocorrer erros na represen¬ 
tação do fenômeno da natureza que estivermos analisando. 

Os problemas ambientais geralmente são complexos e o modelo matemá¬ 
tico é uma aproximação do problema real, representado por expressões mate¬ 
máticas que muitas vezes necessitam de algumas simplificações para obtermos 
uma solução aproximada que resultam em erros. 

1.3 Erros na fase de resolução 

São erros devido ao fato dos equipamentos computacionais terem capacida¬ 
de limitada para armazenar os dígitos significativos de valores numéricos, utilizados 
nas operações elementares de adição, multiplicação, subtração, divisão, etc. 

1.4 Erros na mudança da base 

Considerando que a maioria dos computadores representa os valores nu¬ 
méricos na base binária quando são armazenados, estes são transformados 
em geral, da base decimal para a base binária ou em outra representação. Essa 
transformação pode ser acometida de erros, devido à limitação da representa¬ 
ção do equipamento computacional que estamos utilizando para o processa¬ 
mento dos dados numéricos. 

Dado um número real N, é sempre possível representá-lo em qualquer 
base b, da seguinte forma: 


N > : = Z a i b ' 


em que a, e {0,1,2,3,...,(b-1 )}, com n e m inteiros. 

Exemplo 1.2 

Números escritos na base 2 

a) (1011) 2 =1*2° + 1 *2 1 + 0*2 2 + 1*2 3 

b) (111.01) 2 = 1 *2 2 +0*2- 1 + 1*2° + U2 1 + 1 *2 2 




Exemplo 1.3 

Números escritos na base decimal 

a) (231 ) 10 =1*10° + 3*10 1 + 2*10 2 

b) (231.35) 10 = 5*10 2 + 3*10' 1 + 1*10° + 3*10 1 + 2*10 2 

Mudança da base decimal para a base binária 

Número na base decimal com somente a parte inteira. 
Procedimento padrão: Divisões sucessivas. 

Exemplo 1.4 

Podemos escrever o número (25) 10 na base 2, como segue: 

(25) 10 = 1*2° + 0*2 1 + 0*2 2 + 1*2 3 + 1*2 4 
= ( 11001) 2 

isto é, aplicando o processo de divisões sucessivas temos: 

25+ 2 = 12 e resto = 1 
12+2=6 e resto = 1 
12+2 = 6 e resto = 0 
6+2 = 3 e resto = 0 
3 + 2 = 1 e resto = 1 


Podemos esquematizar o processo de divisões sucessivas, usando o Método 
da Chave, da seguinte forma: 


25 

1 


12 

0 


2 

6 

0 


2 

3 

1 


2 

1 

1 


2 

0 


Para obter o número 25 na base 2, basta tomar os restos das divisões su¬ 
cessivas, de 25 por 2, do quociente 12 por 2, e assim por diante, do final para o 
início no esquema anterior. 


Assim, escrevemos o número 25 como segue: 


25 = 1*2° + 0* 2 1 + 0*2 2 + 1*2 3 + 1*2 4 
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E, portanto temos: 


(25) 10 ={11001) 2 

Definição 1.1 Erro de Arredondamento 

Dizemos que um número x foi arredondado na posição n se todos os dígitos 
de ordem (n + 1) são desprezados da seguinte forma: 

O dígito de ordem n é acrescido de uma unidade se o de ordem (n +1) for 
maior ou igual a 5. Caso contrário, o número x é representado apenas com os 
n dígitos iniciais. 

Exemplo 1.5 

a) O número x = 0.123577, usando a regra de arredondamento anterior, obte¬ 
mos o número x arredondado com 4 casas decimais da seguinte forma: 


x = 0.1236 


b)0 número x = 123343, usando a regra de arredondamento anterior, obte¬ 
mos o número x arredondado com 3 casas decimais da seguinte forma: 


x = 0.123 


Definição 1.2 Erro de Truncamento 

Quando representamos uma função através de uma série infinita, o erro no 
valor de f(x) ao truncarmos a série após um número finito de termos é chamado 
de erro de truncamento. 

Exemplo 1.6 

a)Consideremos a representação de uma função f(x) utilizando a Série 
de Taylor, nas vizinhanças do ponto x, como segue: 

m = m ++ r g >(x) (x g f ) %... + f<" > (x) (x n f ) " +... 


em que, f (n \x) é o valor da n-ésima derivada da função f(x) no ponto x. 

Quando truncamos a série no 3 a termo, isto é, considerando apenas os 
termos até a derivada de ordem 2, na expressão acima, temos um erro cometido 
nessa aproximação, como segue: 









b) Consideremos o desenvolvimento da função f{x) = e x em Série de Taylor, 
nas vizinhanças do ponto x = 0, isto é: 

Como, f{x) = e x , f (1 \x) = e x , f (2 \x)e x ... 

Ainda, f(0) = 1, f (f \0) = 1, f (2) (0) = 1... temos que: 


e x = 1 + x + — + — + ... + — + ... 


2! 3! 


n! 


Considerando apenas os quatro primeiros termos da série, temos: 


X2 X ^ i 

e x = 1 + x + — + — =- (x 3 + 3x 2 + 6x + 6) 
2! 3! 6 


Para x = 2, temos e 2 = 6.33333, que é um valor com erro absoluto bem 
significativo, quando comparado com o valor e 2 = 7.38906 obtido numa calcu¬ 
ladora científica que armazena uma quantidade maior de termos da série. 

Definição 1.3 Erro Absoluto 

Quando substituímos uma solução exata x de um problema por uma solução 
aproximada x', cometemos um erro chamado Erro Absoluto: 



Definição 1.4 Erro Relativo 

Quando consideramos o Erro Absoluto cometido em relação a uma gran¬ 
deza numérica, chamamos de Erro Relativo: 








Exemplo 1.7 


a) Consideremos o valor exato x = 2345.713 e o valor aproximado 
x = 2345.000. 

Nesse caso temos: 

E a = | 2345.713 - 2345.000 \ = 0.713 e, 

E = = 0713 = 0.00030396 

r | x | 2345.713 

b) Consideremos o valor exato x = 1.713 e o valor aproximado x' = 1.000. 


Então, 


E a =\1.713-1.000\ = 0.713 

E r = = °^ 3 - = 0.416229 

r | x | 1.713 

Observe que nos exemplos anteriores a) e b) o erro absoluto é o mesmo, 
embora o erro cometido pela aproximação seja muito mais significativo no exem¬ 
plo b). No exemplo a), o erro relativo é da ordem de 0.03% enquanto no exemplo 
b) é da ordem de 41,6%. 

Exemplo 1.8 

Considere as distâncias entre as cidades: 

a) São Carlos a Campinas -170 km 

b) São Carlos a Ibaté - 20 km 


Ao ser perguntado sobre qual a distância no item a) e b), uma pessoa in¬ 
forma como a seguir: 

a) São Carlos a Campinas -160 km 

b) São Carlos a Ibaté -10 km 








Qual o erro cometido nessas informações? 


E a =\l70 -160 \ = 10 ema) 


E a = \20 -10\ = 10 emb) 



10 

170 


em a) 



10 _ 

20 


1 

— emb 
2 


Podemos concluir que o Erro Absoluto em a) ou b) é o mesmo, não for¬ 
necendo informações significativas sobre o quanto estamos errando naquela 
informação. 

O Erro Relativo nos informa que ao errar 10 em 170 não foi cometido um erro 
muito grande na informação, pois a grandeza numérica 10/170 = 0.0588 = 5.88%, 
mas ao errar 10 em 20, foi cometido um erro de 10/20 = 0.5 = 50% na informação. 

Podemos concluir que o Erro Relativo nos fornece mais informações sobre 
a qualidade do erro que estamos cometendo num determinado cálculo, uma vez 
que no Erro Absoluto não é levada em consideração a ordem de grandeza do 
valor calculado, enquanto que no Erro Relativo essa ordem é contemplada. 

1.5 Erro Absoluto e Erro Relativo nos procedimentos numéricos 

Como utilizar as definições de Erro Absoluto e Erro Relativo? 

Sabemos que nos procedimentos numéricos em geral, quando resol¬ 
vemos um problema, geramos uma sequência de soluções aproximadas 

^0’ ’■■■ ’^k> Xk+1’■■■’ X- 

Se essa sequência for convergente para a solução x, na medida em que k 
cresce, k °°, as soluções aproximadas x k e x k+1 tendem a ficar próximas. 

As definições de E a e E r definidas anteriormente serão usadas para medir 
o quão próximas estão às soluções x k e x k+1 e dessa forma podemos inter¬ 
romper a sequência x k gerada, utilizando um Critério de Parada, com o Erro 
Absoluto ou o Erro Relativo, como segue: 







Dada uma tolerância por um número £ > 0, suficientemente pequeno, con¬ 
sideramos que: 


E a = \ x k+1 -x k |< e 


E = ' Xk+1 * k I < £ 
x k+1 


Caso E a ou E r seja menor do que o número e, entendemos que as soluções 
x k e x k+1 estão próximas e podemos interromper a sequência tomando o último 
valor calculado como a solução aproximada do problema. 

Exemplo 1.9 


Podemos resolver a equação x -2 = 0 utilizando o processo iterati- 
1 _ 

2 


vo x k+1 = - 


f 2' 

x k + — 


V 


,k = 1,2,.... A medida que variamos k = 1,2,... gera- 


y 


mos uma sequência de soluções aproximadas x 0 ,x v ... que converge para 
x = ^2 = 1.414213562. 

Para uma tolerância dada e = 0.0001 e x 0 =1, uma solução inicial, temos 
a sequência de soluções aproximadas convergente para a solução x = 42, ge¬ 
rada pelo processo iterativo, isto é: 

Tomando x 0 = 1, uma solução inicial qualquer e, no processo iterativo dado, 
variando k = 1 , 2 ,... temos: 


Para k = 0^>x 1 = 1.50000 




= 0.33333 > e 


Parak = 1^x 2 =1.41667 


\x 2 -x , 


= 0.05882 >£ 


Parak = 2->x 3 =1.41422 -> ^ * 2 ' =0.00173 >e 


Parak = 3->x 4 = 1.41421 -> X . 3 ' =0.00001 <e 


'4 


22 











Dessa forma, usando o Erro Relativo como o Critério de Parada, e uma 
vez que este está satisfeito na sequência quando k = 3, temos que a solução 
x = x 4 =1.41421. 

Nota: 

Podemos ainda considerar os erros de representação que, em geral, são 
utilizados nos computadores para representação dos números, chamados de 
ponto flutuante. 

A união de todos os números na representação em ponto flutuante é cha¬ 
mado de Sistema de Ponto Flutuante. Detalhes sobre este tópico podem ser 
vistos na Referência [1]. 

1.6 Propagação dos erros 

A propagação dos erros pode ser observada quando utilizamos um pro¬ 
cesso numérico para buscar a solução de um determinado problema. Esse 
processamento envolve um número muito grande de operações elementares, e 
os erros acumulados no final dos processos numéricos interferem na qualidade 
da solução aproximada do problema. 

1.7 Mapa Conceituai 

Um Mapa Conceituai é um instrumento utilizado no processo de ensino, 
aprendizagem e avaliação, por meio de uma representação gráfica de um con¬ 
junto de significados conceituais, o qual pode vir acompanhado de um pequeno 
texto explicativo e de símbolos que ajudam a organizar, sistematizar, estudar e 
detectar as ideias principais do tópico abordado. 

Por exemplo, apresentamos a seguir um Mapa Conceituai preliminar sobre 
os erros em processos numéricos. 





I 

t 

f 1 

Solução Aproximada 

V_'_/ 

Figura 1.2 Um Mapa Conceituai sobre erros. 
































1.8 Exercícios 


1. Converta os números decimais para a base binária: 
a) 45 b) 2978 c) 7699 


2.Resolver a equação x 2 -7 = 0 usando o seguinte processo iterativo: 


X k+1 ~ ' 


1 


x k +- 


, k = 0,1,... com e =10 4 no critério de parada partindo de 


um ponto inicial x 0 = 1. 

3. Usando a técnica de arredondamento, representar os números abaixo 
com 4 casas decimais: 


a) 0.4567897 b) 67.456789 c) 12.677777 


4. Sabemos hoje que a distância média da Terra à Lua é de 384000 km. A 
distância da Terra à Lua obtida por cientistas da antiguidade, na segunda meta¬ 
de do séc. II a.C., foi de 402500 km. 


a) Calcule o erro absoluto cometido nessas informações. 

b) Calcule o erro relativo cometido. Qual é o valor do erro relativo em 
porcentagem? 

5. Faça um mapa conceituai detalhado sobre os erros em processos nu¬ 
méricos, considerando as definições, introduzindo labeis indicando se aprendeu 
(A), não aprendeu (N), se gostou (G), detestou (D ), achou interessante (I), etc. 

Observação: Os mapas conceituais podem ser feitos com os softwares dis¬ 
poníveis como: Xmind <http://www.xmind.net>, CMAPS <http://cmap.ihmc.us> 
ou mesmo a barra de desenhos do Office, Word ou Powerpoint. 







UNIDADE 2 


Solução Numérica de Sistemas 
de Equações Lineares 



2.1 Introdução 


O conhecimento de métodos numéricos para resolução de sistemas de 
equações lineares torna-se necessário devido ao fato de que muitos problemas 
após sua modelagem matemática correspondente resultam na resolução de um 
sistema de equações lineares. 

Esses sistemas lineares em geral são de grande porte, a matriz A dos coe¬ 
ficientes possui muitas linhas e colunas, e necessitam de métodos numéricos 
estáveis, isto é, métodos em que os erros de arredondamento não destroem a 
qualidade da solução aproximada obtida, e com pouco tempo de processamen¬ 
to nos computadores. 

Os sistemas de equações lineares podem ser quadrados, quando o nú¬ 
mero de linhas é igual ao número de colunas, ou retangulares, caso contrário. 
Denominamos sistemas consistentes quando apresentam pelo menos uma so¬ 
lução, e inconsistentes, caso contrário. 

Para ilustrar os sistemas de equações lineares, considere o problema de 
composição de areias para filtros de purificação de água, conforme Figura 2.1 
a seguir: 



Figura 2.1 Filtros de purificação de água. 

Devido a sua propriedade de grande permeabilidade, a areia é usada na cons¬ 
tituição de filtros de Estações de Tratamento de Águas de abastecimento (ETA) 
como meio filtrante, por interceptar as impurezas existentes na água afluente. 

Para construção das unidades de filtração de uma ETA, uma usina de agre¬ 
gados foi contratada para fornecer a areia destinada à filtração. Para atender ao 
pedido, a usina dispõe de areias brutas provenientes de três portos passíveis de 
exploração com composições granulométricas distintas, conforme Tabela 2.1: 









































Tabela 2.1 


Faixa granulométrica (mm) 

Porto 1 

Porto 2 

Porto 3 

0.42-0.59 

0.28 

0.10 

0.15 

0.59-0.71 

0.12 

0.22 

0.10 

0.71-0.84 

0.14 

0.12 

0.18 


Na Tabela 2.1, tomando 1m 3 da areia bruta proveniente do Porto 1 têm-se 
0.28m 3 de grãos na faixa 0.42-0.59mm. As areias brutas são peneiradas para 
obter as areias nas faixas especificadas. 

Essas areias serão dispostas em camadas que deverão obedecer às com¬ 
posições granulométricas estabelecidas por norma, nas quantidades a seguir, 
conforme Tabela 2.2: 


Tabela 2.2 


Faixa granulométrica (mm) 

Volume de areia (m 3 ) 

0.42-0.59 

20 

0.59-0.71 

10 

0.71-0.84 

12 


Deseja-se determinar quais as quantidades de areia de cada porto devem 
ser extraídas de forma que atendam a demanda. 

Na modelagem matemática do problema, sejam x y , j = 1,2,3 as quantida¬ 
des de areia em m 3 a serem retiradas dos respectivos portos P t , i = 1,2,3. 

Dessa forma, para obtermos as quantidades de areia de cada porto que 
devem ser extraídas de forma que atendam a demanda, temos que encontrar a 
solução do seguinte sistema de equações lineares: 

0.28x +0.10x +0.15x =20 

1 2 3 

< 0.12x +0.22x + 0.1 Ox =10 

1 2 3 

0.14x +0.12x +0.18x =12 

lí 2 3 

Observações 

a) Nesse problema algumas simplificações foram feitas, por exemplo, não 
foram considerados os custos de extração de areia nos portos, e ainda 
consideramos apenas três portos de extração. À medida que aumentamos 
o número de portos, os sistemas terão mais equações e mais variáveis. 

b) A resolução deste sistema de equações, requer o conhecimento de mé¬ 
todos adequados e eficientes de resolução, quais serão apresentados 
nesta unidade, como segue: 

















Considere o sistema de equações lineares Ax = b, em que A=(a jj ), 
i,j = 1,...,n é matriz dos coeficientes, x = (x 1 ,x 2 ,...,x n ) t o vetor incógnita 
transposto, e b = (b p b 2 ,...,b n Y é o vetor dos termos independentes transposto. 
Consideramos que o determinante de A, det( A) * 0, o que nos garante a unicidade 
da solução para o sistema linear. 

Representamos o sistema de equações lineares da seguinte forma: 


a x + a x + ... + a x = b 

11 1 12 2 In n 1 

a x + a x + ... + a x = b 

12 1 22 2 2n n 2 


a x + a x +... + a x 

nl 1 n2 2 nn n 


b 

n 


Na forma matricial temos: 


(a a ... a 


fx ) 


f b } 

11 12 In 


1 


i 

a a ... a 


X 


b 

21 22 2n 


2 

= 

2 

a a ... a 


X 


b 

^ nl n2 nn ^ 


K n ) 


k n y 


A x b 


Resolver o sistema anterior dado consiste em determinar um vetor 
x = (x u x 2 ,...,x n ) t que satisfaça todas as equações simultaneamente. 

Podemos representar graficamente a solução de um sistema no 3i 2 , como 
no exemplo a seguir: 

Exemplo 2.1 

x + x = 2 

1 2 

< 

2x + x =3 
l 1 2 

Observe que det( A) = -1*0 e que a solução x = (1,1) 'encontra-se na in- 
tersecção das duas retas, representando as duas equações lineares, conforme 
Figura 2.2: 












Graficamente temos: 



Figura 2.2 Solução de um sistema linear no plano. 

Podemos observar na Figura 2.2 que a solução do sistema de equações 
encontra-se na intersecção das retas, isto é, no ponto x 

Podemos ainda representar a solução de um sistema de equações no es¬ 
paço 9l 3 , isto é, um sistema de equações com 3 equações e 3 incógnitas, con¬ 
forme a seguir: 

Exemplo 2.2 


2x 

+ 

X 

+ Ox = 

3 

i 


2 

3 


-4x 

+ 

2x 

+ X = 

-2 

i 


2 

3 


Ox 

i 

+ 

4x 

2 

+ X = 

3 

4 


Nesse caso, a solução x = (1 , 1,0) 1 do sistema de equações encontra-se 
na intersecção dos três planos representados pelas equações do sistema dado, 
conforme ilustrado na Figura 2.3 a seguir: 
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Figura 2.3 Solução de um sistema linear no espaço. 








Métodos de Resolução 

Na resolução de sistemas de equações lineares, podemos utilizar métodos 
diretos ou métodos iterativos. 

2.2 Métodos diretos 

Um método direto para resolver o sistema Ax = b, em que A= (a,/), 
i, j = 1n com det(A)*0, consiste em determinar exatamente o vetor solução 
x = (x 1 ,x 2 ,...,x n ) t , se não fossem os erros provenientes do processamento de 
cálculos envolvidos nos algoritmos. 

2.2.1 Sistema triangular inferior 

Seja Ax = b um sistemade equações lineares, em que A= ( a tj ), i,j = 1,...,n 
é triangular inferior, isto é, os seus coeficientes (a /y ) = 0 sempre que i < j e com 
a n *0 para 




Representamos um sistema triangular inferior por: 


a x 

ii i 


a x + a x 

21 1 22 2 



a x + a x +... + a x = b 

_ nl 1 n2 2 nn n n 


Para construir o algoritmo para resolução do sistema triangular inferior, 
destacamos a linha genérica (/'), isto é: 

a i1 X 1 + 3/2 X 2 + ... + 3n X, = bj 


Considerando a /( *0 e isolando a incógnita x, na equação anterior, temos: 


x ( . = (b; - a i1 x 1 + a i2 x 2 + a i3 x 3 +... + a ff _,x M ) 





ou seja, 


^ij x i 
i=i 


Dessa forma, podemos escrever o algoritmo para resolução de sistema 
triangular inferior como segue: 

Algoritmo 2.1 

Faça: 



Para / = 2,...,n, faça: 

i-1 

1 a ij x i 

j = i 


Exemplo 2.3 

Considere o sistema de equações lineares triangular inferior: 


~6 

0 

0~ 

V 


~6 

1 

5 

0 

x 2 

= 

6 

0 

1 

1 

_ X 3_ 


1 


Usando o Algoritmo 2.1 temos: 



_ b 2 -a 21 x 1 _ 6-1(1) _ 1 
3^2 5 

_ b 3 ~ a 31 x , - a 32 x 2 _ 1 - 0(1) - 1(1) _ Q 

ã 33 



1 
















Portanto, temos a solução do sistema x = ( 1,1,0) 1 . 


2.2.2 Sistema triangular superior 


Seja Ax = b um sistema de equações lineares, em que A={a jj ), /', j = 1,...,n 
é triangular superior, isto é, os seus coeficientes (a íy ) = 0 sempre que i> j e 
com a /; * 0, para 

i = 


Representamos um sistema de equações triangular superior por: 
a n x i + a 12 x 2 + ... + a 1n x n = b, 

a 22 X 2 + ■ ■ ■ + a 2n X n = ^2 

< 

a nn X n ~ ^n 


Para a construção do algoritmo destacamos a linha genérica (/), isto é: 


a ii X i + a ii +1 Xj +1 + a i i+2 X i+2 + a ii+3 X i+3 + ■■■ + a in X n - b i 


Considerando a„ *0 e isolando a incógnita x, na equação anterior, temos: 


n 

I a U X i 

j=i+1 


Dessa forma, podemos escrever o algoritmo para resolução de sistema 
triangular superior da seguinte forma: 

Algoritmo 2.2 

Faça: 



x 


n 







Para / = faça: 
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X a ü x j 

i= i+i 


Exemplo 2.4 

Considere o sistema de equações lineares triangular superior: 


5 

1 

0 ~ 

X 1 


~ 6 ~ 

0 

7 

-1 

x 2 

= 

6 

0 

0 

8 

_ X 3_ 


8 


Usando o Algoritmo 2.2 temos: 



_ b 2 ã 23 X 3 _ 6 ( 1 )( 1 ) _ 

2 ~ ~ -7 ~ 



x _ - a 13 x 3 _6 - ( 1 )( 1 ) - ( 0 )( 1 ) _ ? 

? ã n 5 

Portanto, temos a solução do sistema x = {1,1,1)'. 

Observação 

O Esforço Computacional E c de um Algoritmo é a quantidade de opera¬ 
ções elementares envolvidas durante a aplicação do mesmo. 

No caso da solução de um sistema de equações triangular superior ou inferior, 
0 esforço computacional dos Algoritmos 2.1 e 2.2 é E c = n 2 operações elementa¬ 
res, sendo n 0 número de operações de divisão, ^ 0 número de operações 
de adição (ou subtração) e n ^ n — 0 número de operações de multiplicação. 

No exemplo anterior, temos que 0 esforço computacional E c = n 2 ={3 ) 2 =9, 
em que três operações são de divisão, três operações são de adição ou subtra¬ 
ção e três operações são de multiplicação. 

















Experiências computacionais mostram que o tempo computacional envol¬ 
vido nessas operações é pequeno, tornando os sistemas triangulares bastante 
atrativos, os quais serão usados no método de decomposição posteriormente. 

Definição 2.1 Menores principais 

Seja a matriz A = (a ij ), i,j = 1,...,n. Denominamos menores principais de 
ordem k da matriz A por: 

A k =det(A k ) 


em que A k = {a jj ), i,j = 1,...,k é formada pelas k primeiras linhas e k primeiras 
colunas de matriz A. 

Exemplo 2.5 

Considere a matriz: 


A = 


5 

1 

1 


0 

2 

2 


0 

0 

1 


Cálculo dos menores principais: 


A, 






0 

2 

2 


0 

0 

1 


= 10 


Métodos de Decomposição 

Os métodos de decomposição para resolução dos sistemas de equações 
lineares Ax = b consistem em decompor de forma única a matriz A do sistema 
no produto de duas matrizes triangular inferior e superior respectivamente, e 
resolvermos dois sistemas triangulares inferior e superior respectivamente, os 
quais são facilmente resolvidos como os exibidos anteriormente. 


2.2.3 Método de Decomposição LU 


O Método de Decomposição LU consiste na decomposição única da matriz 
A = (a jj ), i,j = do sistema de equações, no produto de uma matriz L = {í ij ), 
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i,j = 1,...,n triangular inferior, com os elementos da diagonal l- A =7, / = 1,..., n por 
uma matriz U = (u /y ), i,j = 1,n triangular superior. 

Dessa forma, enunciamos 0 seguinte resultado: 

Teorema 2.1 

Seja uma matriz A=(a ij ), i,j = Se os menores principais de A, A,.*0, 
/ = 1,2,..., n-1, então Ase decompõe, de maneira única, no produto de uma ma¬ 
triz triangular inferior L=(f ij ), i,j = com £ u =1, por uma matriz triangular 

superior U={u ij ),i,j = 

n 

Além disso, det (A) = det (U) = Y[ u a 

i = i 

Prova: Referência [1]. 

Processo de Decomposição LU 

Por facilidade de entendimento, vamos decompor uma matriz A = (a, y ), 
i,j = 1,...,3, isto é, n = 3 e, em seguida, generalizamos os resultados para qual¬ 
quer dimensão como segue: 

Considere LU = A: 


' 1 

0 

0 ~ 



U 12 

U 13 



a i2 

a i3 

C- 21 

1 

0 


0 

u 22 

U 23 

= 

ã 21 

ã 22 

a 23 

J 3 I 

£32 

1 


0 

0 

U 33_ 


a 31 

a 32 

a 33_ 


LU A 


Dessa forma, temos uma igualdade envolvendo 0 produto de duas matrizes e 
determinamos as linhas da matriz U e as colunas da matriz L da seguinte forma: 

• I a linha de U 

1 u„ = a„ = a 11 

1 u 12 = a 12 —> u 12 = a 12 
1 u 13 = a 13 —» u 13 = a 13 










I a coluna de L 


( ii - 1 


^ 21 U 11 _ a 21 ^ 21 


^ 31 U 11 ~ a 31 ^ 31 



/Un 

/Uu 


• 2 a linha de U 

-é* 21^12 ^22 — ^22 ^ Ü22 — ^22 ^ 21^12 

^21 a i3 + ^23 = %3 ^23 = ^23 — ^ 21^13 


• 2 a coluna de L 

í 2 2 = 1 

C31 ^12 + í 32 a 22 = a 32 ^32 = ( a ,32 — ^31^1 2 ) ^22 

• 3 a linha de U 

^31^13 + t 32 a 23 + a 33 = a 33 a 33 = a 33 ~ ^ 31^13 ~ ^32^23 

• 3 a coluna de L 

^33 = ^ 

Podemos generalizar esse procedimento para uma matriz A = {a ij ), 
i,j = 1,...,n e obtemos: 

Matriz triangular U 

i-i 

u n= a vr^L 'ik u kj U = U.,n i<j. 

k=1 




Matriz triangular L 


i-1 

a ij ~ S ^ik U kj 

4 = -—- i,j = 1,—, n / > j. 

u ü 

Assim, podemos escrever o Algoritmo de Decomposição, A = LU, como segue: 

Algoritmo 2.3 

Para m = 1,...,n-1, faça 

Para j = m, m + 1,..., n faça 

m-1 

u m j=a mJ -^( mk u kJ 

k=1 


Para / = m + 1,..., n, faça 


m-1 


a im ^ ik a, 


km 


f . =- 

im 


i =1 

' mm ' 


k=1 


n-1 

^nri ~ a nn — ^ (■ 1 nk C/, 
k=1 


J kn 


= 1 


Exemplo 2.6 

Considere a matriz A dada por: 


A = 


2 

-4 

0 


1 0 
2 1 
4 2 


Como A ; = Ia„ I = 2*0, A 2 = 



a i2 


2 

1 

a2i 

a 22 


-4 

2 


= 8*0, podemos escrever 


A = LU de forma única e seguindo os passos do Algoritmo 2.3: 












I a linha de U 


• I a coluna de L 

Cu = 1 

P 21 = a 21 / u u = -4 / 2 = -2 
£31 = a31 / ^11 = O / 2 = O 

• 2 a linha de U 

u 22 = a 22 - 1 21 u 12 = (2) - (2)(2) = 4 

U 23 = &23 — P 21 U 13 = 0 ) ~ ( — 2)(0) = 1 

• 2 a coluna de L 

P 22 = 1 

i 32 = (a 32 - i 31 u 12 )/u 22 = ((4) - (0)(1)) / 4 = 1 

• 3 a linha de U 

U 33 = ã33 ~ ^31 U13 — ^32 U23 — (2) — mO) — (1)(1) = 1 

• 3 a coluna de L 

C'33 = 1 

Assim, obtemos as matrizes triangular inferior L e superior U 



1 

0 

0 ' 


~2 

1 

0 ~ 

L = 

-2 

1 

0 

u = 

0 

4 

1 


0 

1 

1 


0 

0 

1 








Aplicação na resolução de sistemas de equações lineares 


Seja Ax = b um sistema de equações lineares, em que A = (a /y ), i,j = 1,...,n 
com det(A) * 0. 

Considerando que a matriz A do sistema satisfaz às condições do Teorema 
2.1, temos que A = LU e, portanto, o sistema Ax = b pode ser escrito: 

(LU) x = b 

Nesse caso, denominando Ux = y, a solução do sistema dado consiste na 
resolução de dois sistemas triangulares (inferior) e (superior), respectivamente, 
como segue: 

í Ly =b 
I ux = y 


Para resolver os sistemas de equações triangulares obtidos, basta usar os 
Algoritmos 2.1 e 2.2, respectivamente. 

Exemplo 2.7 

Considere o seguinte sistema de equações lineares: 


2 

1 0 ' 


V 


3 

-4 

0 

2 1 

4 2 


x 2 

_ X 3 _ 

— 

1 

1 

4^ Í\D 

i_ 


Para aplicação do Método de Decomposição LU, verificamos as condições 
de existência e unicidade da matriz 

A = LU , utilizando o Teorema 2.1: 

1. Temos que: 


A ,= 2=2*0 


A 2 - 


2 

-4 


1 

2 


= 8*0 












det(A) = 


2 1 0 
4 2 1 


= 8*0 


0 4 2 

Portanto, a matriz A satisfaz condições do Teorema 2.1 e det(A) * 0 e po¬ 
demos resolver o sistema dado usando o Método de Decomposição LU. 

2. Construção das matrizes L e U: 

Usando o Algoritmo 2.3 temos a decomposição da matriz A = LU da se¬ 
guinte forma: 



1 

0 

0 ~ 


~2 

1 

0 ~ 

L = 

-2 

1 

0 

e U = 

0 

4 

1 


0 

1 

1 


_0 

0 

1 _ 


3. Solução dos sistemas de equações triangulares 
a) Ly = b sistema triangular inferior 


1 

0 

0 ~ 

V 


3~ 

-2 

1 

0 


= 

-2 

0 

1 

1 

_y 3 _ 


4 


Resolvendo o sistema de equações triangular superior obtido, usando o 
Algoritmo 2.1 temos: 

y 1 =3 

< -2y 1 + 1y 2 = 2 y 2 = 4 

y 2 +y 3 = 4 ^y 3 =° 


Assim temos y = (3, 4, 0) { 

b) Ux = y sistema triangular superior 


" 2 

1 

0 


x i 


3 


x i 


1 

0 

4 

1 


x 2 

= 

4 

— > 

x 2 

= 

1 

0 

0 

1 


x 3 


0 


x s 


0 
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Resolvendo o sistema de equações triangular superior, usando o Algoritmo 
2.2 temos: 

x 3 =0 

< 4x 2 + x 3 = 4 — > x 2 = 1 
2x 1 + x 2 = 3 —> x 1 =1 

Portanto, a solução do sistema é dada por x = ( 1,1,0)’. 

Métodos de Eliminação 

Os métodos de eliminação para resolução de um sistema de equações 
lineares Ax = b , consistem em transformá-lo num sistema equivalente triangular 
superior aplicando operações elementares sobre linhas como multiplicar uma 
linha por constante diferente de zero e subtrair de outra linha, as quais serão 
exibidas no método a seguir. 


2.2.4 Método de Eliminação de Gauss 


Considere o sistema de equações lineares Ax = b, em que A = (a jJ ), 
i,j = 1,...,n é matriz dos coeficientes, com det(A) * 0, x = (x p x 2 ,xj é o ve¬ 
tor incógnita transposto e b = (£>,, b 2 ,...,b n ) t é o vetor dos termos independentes 
transposto. 

Representamos o sistema de equações lineares da seguinte forma: 


a„x 7 + a 12 x 2 + ... + a 1n x n = b, 
a 21 x 1 + a 22 x 2 + ... + a 2n x n = b 2 


a m x i + a m x 2 + 


+ a nn X n=b n 


O Método de Eliminação de Gauss, com pivotamento sobre os elementos 
da diagonal da matriz A, consiste em transformar o sistema dado por meio de 
operações elementares sobre as linhas, em um sistema equivalente triangular 
superior, tomando, em cada passo, como pivô, os elementos diferentes de zero 
da diagonal da matriz A. 


(A b) 









em que A (n 1) x = b (n 1] é um sistema triangular superior depois de aplicados 
n -1 passos. 

Consideremos o sistema dado, escrito na seguinte forma: 


, 0 ), 








a i 1 X 1 + a i 2 X 2 + a i 3 x 3 + a 14 x 4 + ... + a 1n x n - a 1n+1 










a 21 X 1 + a 22 X 2 + a 23 X 3 + 3 24 X 4 + ••• + ã 2n X n ~ S 2n+1 


~U) V , _(Í)„ , _(f) v 

d ní A í + d n2 X 2 + «n3 A 3 


, a (í) v _ _(<) 

• • ■ ‘ cl nrv / ' n ct nn+1 


Considere a matriz aumentada, isto é, acoplamos o vetor b dos termos 
independentes juntamente com a matriz A do sistema, para que este também 
sofra simultaneamente todas as operações de pivotamento, de forma a obter no 
final dessas operações um sistema triangular superior equivalente ao sistema 
dado inicialmente. 

Assim, podemos escrever a matriz aumentada da seguinte forma: 



\a (1) 

a ii 

a {1) 

d 12 

a {1) 

d 1n 

a {1) 

d 1n+1 

(Ab) = 

a [1) 

d 21 

a (1) 

d 22 ... 

a {1) 

d 2n 

a U) 

d 2n+1 


a {1) 

L d n1 

a {1) 

d n2 ■ ■ ■ 

a (1) 

d nn 

a {1) 

d nn+1 

em que 

a {1} - 

a ü ~ 

%, i,j = 1,. 

.., n+1 

p a {1) 
u d in+1 


Para facilidade de entendimento no procedimento de aplicações de opera¬ 
ções elementares no sistema de equações, também chamado de processo de 
pivotamento sobre linhas, consideramos um sistema de equações com dimen¬ 
são n = 3 e, finalmente, generalizamos operações para um sistema de equa¬ 
ções de qualquer dimensão. 

Considere a matriz aumentada de um sistema de equações lineares Ax = b 
com dimensão n = 3: 







Assim, podemos escrever: 



\a [1) 

a ii 

a {1) 

d 12 

a (1) 

d 13 

■ a (1) ~ 

d í4 

(At >) = 

a (1) 

d 2i 

a (1) 
d 22 

a (1) 

a 23 

• a (1) 

a 24 


a (1) 

a 31 

a {1) 

d 32 

a (1) 

d 33 

■ a 34_ 

em que 

a [1) - 

a ij - 

-- a ij,Í,J 

= 1,.. 

.,4ea, ( 4 


Passos do processo de eliminação 

• Passo 1 

Nesse passo, considere o elemento pivô * 0. Caso este seja igual a 
zero, troca-se as linhas na matriz aumentada, de forma a obter um elemento 
diferente de zero como pivô. 

Considere os multiplicadores do I a passo: 

a (1) 

a u 


(V) ( 0 \ 

e as operações de pivotamento sobre as linhas de [A, b) do tipo: 


a (,) 

a) Multiplica-se a I a linha fixa da matriz aumentada pela constante —e 

a ií 

subtrai da 2 a linha; faz-se a substituição na 2 a linha por esta modificada 
dessas operações. 


b) Multiplica-se a I a linha pela constante 


j ( í ) 

[31_ 

j(í) 

l n 


e subtrai da 3 a linha; faz-se a 


substituição na 3 a linha por esta modificada dessas operações. 


Dessa forma, conservamos a I a linha e tornamos nulos todos os elementos 

da I a coluna abaixo da diagonal na matriz aumentada (A m , b {1) ), pois ao efetuar 
essas operações de pivotamento temos: 




_(í) _ a 21 
d 21 


a {1) - 0 
( i > a n ~ u 


a d) _ a 31 
d 31 


a {1) - D 
(D a n ~ u 


l ii 


l ii 










Efetuando o mesmo cálculo para toda a linha, obtemos os elementos modi¬ 
ficados, os quais são marcados como ap , e realizamos o mesmo procedimento 
para os elementos das linhas abaixo da diagonal: 

Assim, escrevemos a matriz aumentada modificada como segue: 


(A {1 \b m ) 


p 

1 12 

a {1) ■ 

d 13 

a [1) 

d 14 

» ( 2) 

a (2) - 

a (2) 

‘22 

a 23 

a 24 

,(3) 



l 32 

a 33 

a 34 


De forma geral, para um sistema Ax = b , em que A = (a íy ), i,j = te- 

a (1) 

mos a (2) = ap - m j1 a ( ,J\ m i1 = / = 2,...,n, j = 1, 2,... , n + 1 . 

a u 

• Passo 2 

Considere agora o coeficiente a (2 J * 0 como elemento pivô. Caso ap = 0, 
troca-se com as linhas seguintes até que o coeficiente a% ] que ocupa a segun¬ 
da linha e segunda coluna seja diferente de zero. 

Considere os multiplicadores do 2 2 passo: 


p (2) 

a i2 


e faça as operações sobre as linhas na nova matriz aumentada ( A (1) , b (1) ) modifi¬ 
cada com as operações de pivotamento do I a passo: 

a (2) 

a) Multiplica-se a 2 â linha fixa por e subtrai da 3 a linha; faz-se a subs- 

d 22 

tituição na 3 a linha por esta modificada dessas operações. 

Dessa forma, tornamos nulos todos os elementos da 2 a coluna abaixo da dia¬ 
gonal na matriz (A {2 \ b {2) ), pois ao efetuar essas operações nessa matriz temos: 


a 


( 2 ) 

32 


,( 2 ) 

*32 „(2) 


O ' 

,( 2 ) d 22 

122 


= 0 


Efetuando o mesmo cálculo para toda a linha, temos os elementos modifi¬ 
cados, os quais serão marcados como af\ e procedemos igualmente para as 
linhas abaixo da diagonal: 
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Assim, podemos escrever a matriz modificada da seguinte forma: 



1 

a {1) 

d 12 

a {1) 

d 13 

a (1) 

d 14 

(Á 2 \b (2) ) = 

0 

P (2) 

d 22 

a (2) 

d 23 

a (2) 

d 24 


0 

0 

a {3) 

d 33 

a í3) 

d 34 


De forma geral, para um sistema Ax = b, em que A=(a iJ ), i,j = 1,...,n , 
temos: 


o(3> 

d i j — d i j rn i2 d 2j 


( 2 ) 


m, 


i (2) 

*Í2 

3 l2) 

*22 


3,...,n j = 2 ,...,n + 1 


Assim, depois de executados 2 passos, obtemos o sistema inicial dado 
Ax = b na forma equivalente triangular superior, da seguinte maneira: 




^b 




a i1 X 1 + a i2 X 2 + a i3 X 3 


4 2 b 


-,( 2 b 


- a {1) 

- a 14 

,( 2 ) 


a 22 X 2 + a 23 X 3 — a 24 


,( 3 ), 


,( 3 ) 


a 33 X 3 — a 34 


A solução é dada, conforme anteriormente, usando a expressão geral para 
a solução de um sistema triangular superior: 


X: = 


i = i+i 


a ij x i 


a a 


i = 3, 


1 


Generalizando essas operações para um sistema com n linhas e n colunas, 
necessitamos de n-í passos para transformá-lo em sistema equivalente trian¬ 
gular superior. Essas operações são detalhadas no Algoritmo de Eliminação de 
Gauss, descrito como segue: 









Algoritmo 2.4 


a) Construção do sistema triangular superior equivalente 
Para k = 1 faça 

Para / = k + 1,...,n, faça 


m 


a (k) 
m = _ 

ik oW 
d kk 


Para j = k,..., n + 1 


+ 1) _ q(^) 


ij 


U 


_ (k) (k) 

'"ik d kj 


b) Calcular a solução do sistema triangular superior 

Usar o Algoritmo 2.2 para resolução de sistemas de equações triangular 
superior. 

Exemplo 2.8 

Considere o seguinte sistema de equações lineares: 


"3 

2 

1 

~ x i~ 


"4" 

1 

1 

1 

x 2 

= 

1 

2 

1 

1 

. X 3_ 


2 


Como det (A) = 


3 

1 

2 


2 

1 

1 


1 

1 

1 


= 1*0, o sistema possui uma única solução, 


usando o Método de Eliminação de Gauss, temos: 

Considere a matriz aumentada, conforme a seguir: 


(A b) = 


3 2 1 
1 1 1 
2 1 1 


4 

1 

2 


Executando o Passo 1 e o Passo 2, respectivamente, temos: 

• Passo 1 

Considerando que o pivô afl = 3*0, podemos calcular os multiplicadores: 
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m„ = 


j(7) 

*21 

,(?) 

*11 


1 _ 

3 


e m 


21 


a (1) 

_ d 31 _ 

.(») 

*11 


2 _ 

3 


Multiplica-se a I a linha por 1/3 e subtrai da 2 a linha; faz-se a substituição 
na 2 a linha, em seguida multiplica-se a I a linha por 2/3 e subtrai da 3 a linha e, 
novamente, faz-se a substituição na 3 a linha. 

Assim temos a matriz aumentada modificada da seguinte forma: 


(A { 1 \b {1) ) = 


3 

0 

0 


2 1 
1/3 2/3 

-1/3 1/3 


4 

-1/3 

-2/3 


• Passo 2 

Como o pivô a 22 = 1/3 * 0, calculamos os multiplicadores como segue: 


m 32 ~ 


,( 2 ) 

d 32 


,( 2 ) 


-1/3 

1/3 


= -1 


Multiplica-se a 2 a linha por -1, subtrai da 3 a linha e faz-se a substituição na 
3 a linha. Temos a matriz aumentada modificada como segue: 


(A { 2 \b (2) ) = 


3 

0 

0 


2 1 

1/3 2/3 
0 1 


4 

-1/3 

-1 


Assim, obtemos um sistema de equações triangular superior, o qual pode 
ser reescrito da seguinte forma: 

3x 1 + 2x 2 +1x 3 =4 

1/3x 2 + 2/3x 3 = - 1/3 
1x 3 = -1 


Resolvendo o sistema triangular superior equivalente, usando o Algoritmo 
2.2, temos: 

1x 3 — — 1 —^ x 3 = — 1 

< 1/3x 2 + 2/3x 3 = -1/3 -» x 2 = 1 
3x 1 + 2 x 2 +1x 3 =4 x 1 = 1 














Portanto, a solução do sistema é dada por x = ( 1,1, -1) { . 

Observação 

Para verificar-se que a solução obtida é solução do sistema original dado, 
basta substituí-la nas equações do sistema e verificar a igualdade. 

Aplicação 

Uma aplicação de resolução de sistemas de equações lineares é quando 
calculamos a matriz inversa de uma matriz A, em que A = (a jj ), i,j = 1,... ,n. Nesse 
caso necessitamos da resolução de n sistemas de equações lineares, em que 
usamos os métodos diretos vistos anteriormente. 

2.3 Inversas de matrizes 


Considere A = (a i ,), i,j = 1,... ,n uma matriz não singular, isto é, det(A) * 0. 
Então existe uma única matriz A 1 chamada de inversa de A, de tal modo que 
AA -1 =1, em que a matriz I é uma matriz diagonal, conhecida como a matriz 
identidade, cujos elementos da diagonal são todos iguais ao número 1. 

Dessa forma temos: 


f *\1 

a 12 . . 

• 

( X 

A 11 

X 12 ■ 

■ X 1n 


f 1 

0 

0 . 

. 0 ' 

a 2i 

a 22 ■ ■ 

• a 2n 

X 21 

X 22 ■ 

■ X 2 n 

= 

0 

1 

0 

. 0 

y^nl 

a n2 ■■ 

^nn j 

K X n1 

X n2 •• 

• X nn , 


° 

0 

0 . 

. 1 

J 


Portanto, para determinar as n colunas da matriz inversa A 1 , temos que 
resolver n sistemas de equações lineares, usando qualquer um dos métodos 
diretos vistos anteriormente. 

Para calcular a I a coluna da matriz inversa A' 1 , resolvemos o seguinte sis¬ 
tema de equações lineares: 



a 


12 


a 


21 


a 


22 



a 


n2 


In 


2n 


'X ' 
11 


21 


f 1 ^ 

0 


o 

v y 
















Assim sucessivamente, calculamos a n-ésima coluna da matriz A \ resol¬ 
vendo o seguinte sistema de equações lineares: 


a a 

11 12 


a a 

21 22 



(x 

In 


\ 



'o N 

0 


v 


a 


nl 


a 


n2 


a 

nn 


X 

\ nn J 


y 


Exemplo 2.9 

Considere a matriz A dada por: 


A = 


1 

1 

0 


0 2 
1 0 
0 2 


Temos que a matriz inversa de A pode ser escrita como: 


X 

X 

X 

11 

12 

13 

X 

X 

X 

21 

22 

23 

X 

X 

X 

31 

32 

33. 


Usando o Método de Eliminação de Gauss, podemos calcular a inversa da 
matriz A dada, como segue: 


Como det(Á) 


1 0 2 
1 1 0 
0 0 2 


= 2*0 temos que a matriz inversa existe e é única. 


Considerando que AA 1 =1, podemos escrever: 


1 

0 

2 " 


1 

1 

0 


0 

0 

2 



xxx 

11 12 13 

XXX 

21 22 23 

XXX 

31 32 33. 


1- 

1_ 

0 

1 - 

O 

= 

0 

1 

0 

-1 

1 

O 

0 

1 


Temos três sistemas de equações lineares para serem resolvidos, como 
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seguem: 






















I a coluna da matriz inversa de A, 


1 0 2 


X 


V 

1 1 0 


ii 

X 


0 


21 


0 0 2 


1 

* 

Co 

1_ 


0 


2 a coluna da matriz inversa de A, 


~1 0 2 


X 


~0~ 



12 



1 1 0 


X 

— 

1 



22 



0 0 2 


X 

L 32 J 


0 


3 a coluna da matriz inversa de A, 


1 0 2 


X 


" 0 " 



13 



1 1 0 


X 

— 

0 



23 



0 0 2 


X 

L 33 J 


1 


Construímos inicialmente a matriz aumentada (A, I) formada pela matriz A 
e os termos independentes dos sistemas obtidos anteriormente. Aplicamos as 
operações do Método de Eliminação de Gauss, transformando a matriz dada na 
forma triangular superior, e os vetores dos termos independentes são modifica¬ 
dos, como segue: 


' 1 

0 

2 

1 

0 

0 ~ 


' 1 

0 

2 

1 

0 

0 ~ 

1 

1 

0 

0 

1 

0 

~ 

0 

1 

-2 

-1 

1 

0 

0 

0 

2 

0 

0 

1 


0 

0 

2 

0 

0 

1 


Assim, podemos resolver os sistemas triangulares obtidos: 


' 1 0 2 


X 


1 



ii 



0 1 -2 


X 

— 

-1 



21 



0 0 2 


X 

L 3/J 


0 


Resolvendo o sistema de equações triangular superior obtido, temos a so¬ 
lução da 1 s coluna da matriz inversa de A, dada por: 


(x„ x 21 x 31 )’ = (1,-1,0) 1 
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Para determinar a 2 a coluna da matriz inversa de A, resolvemos o sistema 
de equações triangular superior obtido: 


~1 0 2 


X 


~0~ 



12 



0 1 -2 


X 

— 

1 



22 



0 0 2 


1 

* 

l\3 

1_ 


0 


cuja solução é dada por: 

(Xy 2 X22 X 32 ) = (Qt 0) 

Finalmente para determinar a 3 a coluna da matriz inversa de A, resolvemos 
0 sistema de equações triangular obtido como segue: 


1 0 2 


X 


~ 0 ~ 



13 



0 1 -2 


X 

— 

0 



23 



0 0 2 


X 

L 33 J 


1 


cuja solução é dada por: 

{x 13 x 23 x 33 y = (- 1,1,1/2)' 


Portanto, temos a matriz inversa A _1 da matriz dada por: 


A 1 


1 0 -1 

-1 1 1 

0 0 1/2 


Observação 

Quando usamos 0 Método de Decomposição LU no cálculo da matriz inver¬ 
sa, transformamos a matriz A no produto de duas matrizes triangulares, inferior 
e superior respectivamente. Nesse caso, teremos três sistemas de equações li¬ 
neares para serem resolvidos e, consequentemente, dois sistemas de equações 
triangulares inferiores e superiores em cada um deles. 


















2.4 Métodos iterativos para resolução de sistemas de equações 
lineares 

Podemos resolver um sistema de equações lineares Ax = b utilizando para 
isto métodos iterativos, os quais consistem em determinar uma sequência de 
soluções aproximadas, em que cada solução obtida depende da solução ante¬ 
rior pela aplicação do mesmo procedimento. 

Para construir um método iterativo, escrevemos o sistema de equações 
lineares Ax = b, na forma equivalente x = Hx + g. Essa forma equivalente será 
detalhada no Método iterativo de Jacobi-Richardson, conforme segue. 

De maneira geral, obtida a forma equivalente x= Hx + g, a partir de uma 
solução aproximada inicial x (0) , determinamos a sequência de soluções aproxi¬ 
madas considerando o processo iterativo: 

x (í) =Hx (0] +g 
x (2) = Hx (1) +g 


x (k+1) = Hx {k) +g, k = 0,1,2,... 

Nesse caso, temos H = (õ /y ), i,j = 1,..., n chamada de matriz iterativa e um 
vetor g = (g,), i = 1,...,n. 

Assim, a partir de uma solução aproximada inicial x (0) = {x\°\ xx [ n 0) ) 
para a solução exata x = (x 1 ,x 2 ,...,x n ) do sistema Ax = b, determina-se uma 

sequência de vetores x {1 \ x {2) , x°\... e espera-se que seja convergente para a 
solução x, isto é, 

jim x (k) = x 


Apresentamos a seguir alguns resultados e definições necessários para 
entendimento do Método iterativo de Jacobi, como seguem: 

Definição 2.2 Norma de vetores 

A norma de um vetor xeV (espaço vetorial) é uma função ||. ||; V -»9t + , 
satisfazendo às seguintes condições: 

a) li x II > 0; li x II = 0 <^> x = 0 , V x e V 




b) I ax I = | a II x ; Vae 9Í, Vxe V 

c) ||x + y ||<||x|| + ||y ||;Vx, ye V 


Exemplos de normas: 

Considere o espaço 9í n formado pelos vetores x = (Xp x 2 ,..., x n ): 
Temos as seguintes normas de vetores: 

... . i i n | i 

X II = X + X + ... + X =Y X 
11 / | 1 | | 2 | | n | 4- J \ i | 




( n „ X 

X 1 =lx 2 + X 2 + . 

II 2 \ 1 2 

, + x 2 = 

n 

Xh,l 



V’ 1 7 


X 


máx ix 

, x L 

.. Ix f 

[ = máx x 

11 1 1 

1 2 \ 

1 n [ 

\ 1<i<n i \ 


Essas normas são chamadas de norma 1, norma 2 ou Euclidiana e norma 
infinita, respectivamente. 

Observação 

Podemos interpretar as normas definidas anteriormente no 3i 2 , bastando, 
para isso, representar o ponto x = (x 1 ,x 2 ) t no plano e observar, por exemplo, que 
a norma 2 é a distância entre o ponto x e a origem do sistema de coordenadas. 

Exemplo 2.10 

Considere x = (x p x 2 , x 3 ) = (-1,2,-4)& 9í 3 . 

Calculando as normas do vetor x, obtemos: 


x 



-1 + 2 + 


4 


= 7 



( 3 


2 ' 

l x L = 


X 

/ 



l íW 


) 


4l 2 +2 2 +4 2 




II II 

, r 

X [ 

o = maxl 



} = máx{1,2,4} = 4 


Definição 2.3 Sequência convergente 

Dizemos que uma sequência de vetores x (,) = (x p x 2 ,..., x n ) 0) e 9T con- 




















verge para o vetor x = (x p x 2 ,..., x n ) e Si n se 


x (/) -x 


->0, quando /—>°°, para qualquer norma em 9í n . 


Exemplo 2.11 

Considere a sequência x (,) = 


-,0,0,...,0 
i 


e9t n , i = 1,2,... 


V V 

Temos que, à medida que /—>°°, a sequência converge para o vetor 
x = (0,0,..., 0). 


Ou ainda, 
Como 


x (/) -x 


= -, temos 


x (í) -x 


0, quando / —^ 00 . 


Definição 2.4 Norma de matriz 

Considere o espaço 9í(n, n) das matrizes quadradas A = (a,.) /, / = 7,..., n. 
Definimos as seguintes normas de matrizes: 


I A II =M|L = máxY 

i M7 " no i<j< n “ 


I A II = II /A II, = máxY 


1<i<n 


H 


Essas normas são chamadas norma coluna e norma linha, respectivamente. 

Exemplo 2.12 

Considere a matriz A: 


A = 


7 2 3 
1 1 1 
1 2 4 


Então, 


A 



máx y 

1<i<3 ^ 


a 

ij 


máx\3, 5, 8} = 8 

1<j<3 L J 


























>A|| = II /A II = máxY a = máx\6,3,7} = 7 

Moo M ML 1<i<3 ' " islsn L J 


H 


1<i<3 


Definição 2.5 Normas consistentes 

Considere uma norma de um vetor xe9T e uma norma de uma matriz 
Ae 9í(n, n). Dizemos que essas normas são consistentes se: 

I Ax I < I A II x ||,Wte 9?(n,n) e Vxe9í n . 

Propriedade 

||/Ae||<||/A||||e||, VASe 9í(n,n). 

Definição 2.6 Matriz estritamente diagonalmente dominante 

Dizemos que uma matriz A = {a jj ), i,j = 1,...,n é estritamente diagonal¬ 
mente dominante se: 



n 


a 

// 

V 

M 

a 

ij 

H. 


'W 


Exemplo 2.13 


A = 


-4 

1 

1 


1 

-2 

1 


0 

0 

3 


Usando a definição 2.6, temos que: 

/ = /-»|-4|>|í| + |0| 

/ = 2 -> | - 2 | > | /1 +1 0 | 

/ = 3—»|3|>]f| + |f[ 

Portanto, concluímos que a matriz dada é estritamente diagonalmente 
dominante. 












2.4.1 Método iterativo de Jacobi-Richardson 


Considere o sistema de equações lineares Ax = b, em que A = (a /y ), 
i,j = é a matriz dos coeficientes com det(A) * 0 e com os elementos da 
diagonal principal a„ * 0, i = 1,...,n, x = {x 1 ,x 2 ,...,xj é o vetor incógnita trans¬ 
posto e b = (b p b 2 ,... ,b n )’ é o vetor dos termos independentes transposto. 

Representamos o sistema dado na seguinte forma: 

a n x i + a i 2 x 2 + a 13 x 3 + ... + a 1n x n = b, 
a 2 i x i + a 22 x 2 + a 23 x 3 + ■ • ■ + a 2n x n = b 2 

< 

a n1 X 1 + a n2 x 2 + a n3 X 3 + - + a nn X n = b n 

O Método iterativo de Jacobi consiste em escrevermos o sistema Ax = b 
na forma equivalente x = Hx + g, como segue: 

Dividindo cada linha pelo elemento da diagonal e explicitando x 1 na 1 a equa¬ 
ção, x 2 na 2 a equação, x 3 na 3 a equação e até x n na n-ésima equação, temos: 

X 1 = — ( b 1 - a i2 X 2 ~ a i3 X 3 ~ ■■■ ~ a in X n ) 

a i1 

X 2 = ( b 2 — a 21 X 1 — a 23 X 3 — ■ ■ ■ — a 2n X n ) 

a 22 

x n = ^~ (b n -a n1 x 1 - a n2 x 2 - ...-a nn _ f x„_,) 

l a nn 

Na forma matricial temos: 


V 


f 

0 

a 21 

a 31 

_ a n1 

\ 

V 


'A. 1 



a i1 

a i1 

a i1 




a n 

X 2 

= 

a 21 

a 22 

0 

a 23 

a 22 

a 2n 

a 22 


X 2 

+ 

È2_ 

ã 22 



_ a n1 

a n2 


0 




bn_ 



v a nn 

a nn 


/ 



y a nn À 














Assim, podemos escrever x = Hx + g , em que 


h a=< 


0, i = j 


a ij ■ ■ 

i*J 

a,-,- 


é chamada de matriz iterativa e g j = —, / = 

3// 


Podemos escrever o processo iterativo de Jacobi como segue 
x (k+1) = Hx ik) +g k = 1,2,... 

Assim: 


íxS k+1 '< \ 


x 2 ( k+1 ) 


x (k+í) 
v n y 


0 

a 27 

a 31 

_ a n1 

a 7í 

a 77 

a i 1 

a 27 

0 

a 23 

a 2n 

a 22 

a 22 

a 22 

fnl 

ã n2 


.. 0 

a nn 

a nn 



/ K A 


'x/*) A 


x 2 (' < ) 


X„(*> 

v n y 


bi 


*22 


v a ""y 


Podemos, ainda, escrevê-lo na seguinte forma: 

x/* +í ) = 


a 7? a í/ 


. - ^x „(*) + A_ 


-*7 7 


-^7 7 


x 2 ( - k+1 '> = 


- - ^X 3 (*) 

a 22 a 22 


X (k+7) = - ^-x,( k ) - — XJ k 1 - 


hn x (k) + A_ 

n 

a 22 a 22 


- ^= L x n /*)+ A_ 

a n/7 " a nn 















Estudo da convergência 

Teorema 2.2 

Seja uma norma matricial consistente com alguma norma vetorial e 
x (0) e9T uma solução aproximada inicial qualquer para o sistema de equa¬ 
ções. Se I H \\<1, então a sequência de soluções aproximadas definida pelo 
processo iterativo x [k+1) = Hx {k) +g, k = 0, 1 , 2 ,... converge para a solução x 
do sistema de equações Ax = b. 

Prova: Referência [1]. 

Resultado: 

É possível provar que se a matriz A do sistema Ax = b for estritamente 
diagonalmente dominante, então o Método de Jacobi-Richardson gera uma se¬ 
quência de soluções aproximadas convergente para a solução do sistema, pois 
sempre temos || H < 1. 

Apresentamos de forma detalhada o Algoritmo do Método de Jacobi-Ri¬ 
chardson, observando que o critério de parada na sequência gerada de vetores 
é o erro relativo visto anteriormente, apenas usamos a notação de norma, pois 
estamos trabalhando com sequência de vetores. 

Algoritmo 2.5 

a) Forneça uma solução inicial aproximada x (0) = (x\ 0) , x { 2 0) ,..., x ( n 0) ) e e > 0 
uma tolerância fixa. 

Faça k = 0 e PARE = FALSO 

b) Construção da sequência de soluções aproximadas: 

Enquanto PARE = FALSO, faça: 

Para / = /,..., n faça: 


x 


(k+1) 


f í ~ a ij X T + I - a a x 

i=i j=i+i 


\ 


( k ) 
i 


J 


/a ii + J- 
a ii 


X (k+1 ) _ x(k) 

Se -—TT-TT— ü-<8, então PARE = VERDADE 

X (k+1) 


Senão k = k + 1 
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Exemplo 2.14 

Considere o seguinte sistema de equações lineares: 


~4 

1 

1~ 

x i 


" 3" 

1 

4 

1 

x 2 

= 

0 

_1 

1 

5_ 

_ X 3 _ 


-4 


Usando o Método iterativo de Jacobi, podemos a partir de uma solução 
aproximada inicial x (0) = (x y (0) , x { °\ x ( 3 0) )' =(0,0,0)' calcular uma solução apro¬ 
ximada para o sistema dado com uma precisão e = 10'. 

Assim, a construção da matriz iterativa H é dada por: 


H = 


0 

1/4 

1/5 


-1/4 

0 

-1/5 


-1/4 

-1/4 

0 


e para verificar a condição de convergência temos: 


H II = II H II = máx Y I h:, I = máx{2/4, 2/4, 2/5} = 0.5<1 

" L 11 No ° / 1 1 / L J 


Portanto, a sequência de soluções aproximadas é convergente para a so¬ 
lução do sistema dado. 

Cálculo das iterações: 


^(k+1) _ 




y(k + 1) _ 

- 


X 


(k+1) _ 
3 


X 


( k) 


— X 


( k) 


— X 


( k ) 


— X 

4 


( k) 


— X 

4 


(.k) 


— X 


( k ) 


3 

+ — 

4 

+ 0 

4 

5 


k = 0,1,2,... 


Para k = 0 e tomando x (0) = (0,0,0 )', temos x 


|x (í >-x<°> 
II yW II 


(?) _ 


V 


-, 0 ,-~ 

4 5 
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Teste de parada: 


= 1>10 
















Para k = 1, temos: 


x {2) = 0.9500 
xf = 0.0125 
xf = -0.9500 

Portanto, 

x< 2 > = ( 0.9500, 0.0125,-0.9500Y 

II x < 2 ) -x {1) II 

n-n—^ = 0.2000 > 10- 1 

x< 2 > 


x\ 2) = 


A2) _ 


A2) - 


1 (í) 1 ( 1 ) 3 

-x‘ 7) -xV'+ — 

4 2 4 3 4 

--x f (,) - -x<’> + o 

4 1 4 3 

- —x {1) - -x (1) - — 

5 1 5 2 5 


—> 


—> 




Assim, sucessivamente, calculamos: 
Para k = 3, temos: 


x ? (3) = 


y(3) _ 
A 2 


y (3)_ 

A 3 


~-x (2) ~-X (2) +- 

4 2 4 3 4 

--x\ 2) --x {2) +0 

4 1 4 3 

-L X (2)_L X (2)_Í 

5 1 5 2 5 


—> 


—> 


—> 


x ? (3) = 0.9844 
xf= 0.0000 
xf =-0.9925 


Portanto, 

x (3) =(0.9844, 0.0000,-0.9925)’ 

I X (3) _x (2) I 

-—TT-tt —— = 0.0400 < 10- 1 , então pare. 

x (3) 


Logo, x = x (3) = (0.9844,0.0000, -0.9925)* é a solução aproximada para 
sistema dado com a precisão e = 10 








Observação 

Podemos observar que a sequência de soluções aproximadas 
x (1) , x (2) , x (3) , x (4 \... converge para a solução do sistema de equações 
x=(1,0,-1) t . 

Observação 

Os métodos iterativos devem ser aplicados na resolução de um sistema de 
equações lineares Ax = b quando este apresentar a matriz A de grande porte e es¬ 
parsa, isto é, uma matriz com muitos elementos em que a maioria deles são nulos. 

Nesse caso, quando determinamos a matriz iterativa H, esta contém muitos 
elementos nulos, o que diminui o esforço computacional. Caso contrário, deve¬ 
mos utilizar os métodos diretos vistos anteriormente. 


2.5 Exercícios 

1. Considere o seguinte sistema de equações lineares: 


'-5 

1 

1 " 

V 


2 

0 

1 

0 

x 2 

= 

1 

2 

3 

2 

.*3. 


5 


a) Verifique as condições de aplicabilidade do Método de Decomposição 
LU para o sistema dado. 

b) Caso possível, resolvê-lo usando o Método de Decomposição LU. 

2. Determinar a inversa da matriz A do exercício 1), usando o Método de 
Decomposição LU. 

3. Resolver o seguinte sistema de equações lineares usando o Método de 
Eliminação de Gauss. 


~0 

6 

-1 

2 

V 


f 

0 

5 

-1 

1 

x 2 


0 

2 

1 

-1 

-1 

*3 


-2 

0 

-1 

-2 

1 

_*4. 


-1 


4. Desenvolva o Método de Gauss-Jordan para resolução de sistemas de 
equações lineares e faça um exemplo. Referência [1]. 


64 
















5. Considere o sistema de equações lineares: 


~0 

1 

1~ 

V 


~0~ 

2 

3 

0 

x 2 

= 

3 

3 

1 

1 

_ X 3_ 


0 


a) Faça trocas de linhas e, caso possível, resolva o sistema usando o Mé¬ 
todo de Decomposição LU. 

b) Determinar a inversa de A, usando o Método de Gauss-Jordan. 

c) Faça trocas de linhas e, em caso de convergência garantida, resolva o 
sistema usando o Método Iterativo de Jacobi com uma solução inicial 
dada e e = 0.1. 


6. Considere o seguinte sistema de equações lineares: 


-9 

2 

2 

V 


"-5" 

1 

3 

0 

x 2 

= 

4 

-1 

1 

3 

X 3. 


3 


Caso haja convergência garantida para o Método Iterativo de Jacobi, resol¬ 
va o sistema dado a partir da solução inicial x(°) = (0 00) e e = 0.1. 

7. Mostre que || x || m < || x || < n|| x ^, x £ R n . 

8. Considere duas cidades C 1 e C 2 com produções de lixo de 8 e 5 tonela¬ 
das, respectivamente. Numa área apropriada existem dois aterros A 1 e A 2 com 
capacidades de processar 7 e 6 toneladas de lixo, respectivamente. Deseja-se 
saber quais as quantidades de lixo devem ser enviadas das cidades C 1 e C 2 para 
os aterros A 1 e A 2 de forma que todo o lixo seja transportado. 

a) Escreva a formulação matemática correspondente do problema dado. 

b) Resolvê-lo usando o Método de Eliminação de Gauss. 

c) Interprete os resultados obtidos. 

d) Resolva o problema dado, com o Software Matlab, usando o Método de 
Gauss-Jordan. 

e) Resolva o problema dado, com o Software Numérico, usando o Método 
de Eliminação de Gauss - Referência [1]. 
















9. Faça um mapa conceituai detalhado sobre solução numérica de sistemas 
de equações lineares, considerando os métodos vistos, introduzindo labeis indi¬ 
cando se aprendeu (A), não aprendeu (N), se gostou (G), detestou (D), achou 
interessante (I), etc. 




UNIDADE 3 


Solução Numérica de Equações 



3.1 Introdução 


Nas diversas áreas científicas, frequentemente deparamo-nos com problemas 
reais envolvendo a resolução de equações, isto é, dada a equação f(x) = 0, dese¬ 
jamos determinar a solução (ou soluções) real x , tal que f(x) = 0. Por exemplo, 
consideremos a equação f(x) = sen(x) - ln(x) + 2 = 0, desejamos determinar a 
solução real x, tal que f(x) = sen(x)-ln(x) + 2 = 0. 

Métodos iterativos são desenvolvidos para determinar aproximadamente 
essa solução real x, embora tenhamos métodos iterativos específicos para de¬ 
terminar a solução x, quando esta é um número complexo. 

Métodos iterativos são apresentados para determinar a solução x, quan¬ 
do esta é um valor real, e para isso, necessitamos de uma solução aproximada 
inicial x 0 , a partir da qual geramos uma sequência {x n }, n = 1,2,... de soluções 
aproximadas, que sob determinadas condições teóricas converge para a solu¬ 
ção desejada x. 

A solução inicial x 0 pode ser obtida através de recursos gráficos, no qual 
localizamos uma vizinhança ou um intervalo [a,b] em que encontra-se a solu¬ 
ção x, conforme exibimos na Figura 3.1: 



Observando a Figura 3.1, vemos que a solução x tal que f(x) = 0 encon¬ 
tra-se onde a função f(x) corta o eixo das abscissas, isto é, um intervalo em 
que os valores da função possui sinais opostos. 

Podemos, então, tomar uma solução inicial x 0 nas vizinhanças dessa raiz, 
isto é, no intervalo [a, b\ para inicializar a sequência de soluções aproximadas du¬ 
rante aplicação dos métodos iterativos que serão apresentados posteriormente. 









Para ilustrar a solução numérica de equações, considere o problema de 
nível de oxigênio em um rio, conforme Figura 3.2 



Figura 3.2 Rio contaminado. 

Sabendo-se que para calcular o nível de oxigênio N (mg/L) em um rio a 
jusante de uma carga de esgoto, podemos usar a seguinte função: 

N(x) = 10-20{e-°- 15x -e-° 5x ) 

em que x é a distância a jusante em quilômetros, conforme Figura 3.3: 



Figura 3.3 Nível de oxigênio de um rio. 

a) Deseja-se determinar a jusante em que o nível do oxigênio cai pela primei¬ 
ra vez até uma leitura de 5 mg/L, observando-se quais níveis de oxigênio 
abaixo de 5 mg/L são prejudiciais aos peixes e ao ecossistema aquático. 

b) Deseja-se determinar a distância a jusante na qual o oxigênio está no 
valor mínimo. Qual a concentração de oxigênio nessa posição? 









Definição 3.1 

Dizemos que x é uma raiz, ou um zero, da função f(x) se f(x) = 0. 

Exemplo 3.1 

Seja f(x) = x 2 - 5 = 0. Temos que as raízes da equação x 2 - 5 = 0 são: 
x = ± yfíT = 2.236067977 , e nesse caso f(x) = 0. 

3.2 Localização das raízes: Método Gráfico 

Para localizar uma vizinhança para a raiz de f(x), traçamos o gráfico de f(x); 
assim, os pontos em que este corta o eixo das abscissas nos informam a raiz (as 
raízes) de f(x). 

Considere o exemplo anterior, f(x) = x 2 -5, e temos as raízes x = ±V5, 
conforme gráfico ilustrado na Figura 3.4: 



Figura 3.4 Raízes de uma função. 

Observando a Figura 3.4, vemos que o gráfico de f(x) permite identificar 
onde estão aproximadamente as raízes de f(x). Nesse caso, temos as raízes 
x = ± 45 = 2.236067977 com f(x) = 0. 

Podemos ainda transformar a equação f(x) = 0 na forma equivalente, na 
igualdade de outras funções /j(x) = f 2 (x). Os pontos de interseção dos gráficos 
/j(x) e f 2 (x) serão as raízes procuradas, conforme exemplo a seguir: 


Exemplo 3.2 

Considere a equação f(x) = 4x-e x =0. 
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Podemos escrever a equação dada na forma equivalente 4x = e x , isto é, 
f 1 (x) = f 2 (x), com f 1 (x) = 4x e f 2 (x) = e x , conforme gráfico da Figura 3.5: 



Figura 3.5 Intersecção de gráficos de funções. 

Observando a Figura 3.5, vemos que a raiz x encontra-se na intersecção 
dos gráficos de /j(x) e f 2 (x). 

3.3 Métodos numéricos para resolução de equações 

Nesta unidade apresentaremos alguns dos principais métodos para resol¬ 
ver numericamente uma equação. 

3.3.1 Método da Bisseção 


Considere a função f: 9? uma função contínua. Desejamos resolver a 

equação f(x) = 0, isto é, determinar uma solução x real tal que f(x) = 0. 

O Método da Bisseção é baseado no Teorema do Valor Intermediário, o qual 
afirma que se uma função é contínua no intervalo [a,b], e satisfaz a condição 
f(a)f(b)<0, valores de f(a) e f(b) com sinais opostos, então existe xe (, a,b ) 
tal que f(x) = 0, isto é, existe ao menos uma raiz no intervalo \a,b], conforme 
ilustrado na Figura 3.6. 







Figura 3.6 Intervalo contendo uma raiz da função f(x). 

O Método da Bisseção consiste em localizar inicialmente um intervalo 
[a,b], em que se encontra a raiz x, e determinar uma sequência de intervalos 
[a/ = 0,1,..., em que r 0 =a e s 0 = b, de forma que a amplitude do intervalo 
numa iteração seja igual à metade da amplitude do intervalo anterior e que sem¬ 
pre contenha a raiz x. 

A sequência de intervalos será calculada até que a amplitude do intervalo 
seja menor do que uma tolerância s preestabelecida. 

As sequências r jt s, e x, são construídas a partir de um intervalo inicial 
[r 0 , s 0 ] tal que f(r 0 )f(s 0 ) < 0, conforme Figura 3.7: 



Figura 3.7 Intervalo inicial contendo uma raiz da função f(x). 














Determina-se o ponto médio do intervalo [r 0 , s 0 ], dado por: 


x 0 =(r 0 +s 0 )/2 

A partir da solução aproximada obtida x 0 , fazemos a seguinte verificação: 
Se f(x 0 ) = 0, então x 0 é uma raiz de f(x) 

Se f{r 0 )f(x 0 ) < 0, então r, = r 0 e s 1 = x 0 
Se f(r 0 )f(x 0 ) > 0, então ^ = x 0 e s 1 = s 0 

Neste momento, temos o novo intervalo \r p s ; ], cuja amplitude é igual à 
metade da amplitude do intervalo [r 0 , s 0 ] e que contém a raiz desejada x. 

O procedimento é repetido novamente, isto é, calcula-se o ponto médio do 
intervalo [r p s r ] dado por x 1 = (r 7 +s 7 ) /2, conforme observado na Figura 3.8. 



Figura 3.8 Intervalo contendo o ponto médio. 

A partir da solução aproximada obtida x p fazemos a seguinte verificação: 
Se f(x 1 ) = 0 , então x 1 é uma raiz de f(x) 

Se f(r 1 )f(x 1 )<0, então r 2 =r 1 e s 2 =x 1 
Se f(r 1 )f(x 1 )>0, então r 2 =x 1 e s 2 =s 1 

Assim sucessivamente, temos um intervalo genérico [r jt s ( ], em que calcu- 
la-se o ponto médio deste intervalo, dado por x, = (ç + s,) /2\ 

Se f(Xj) = 0, então x ( . é uma raiz de f(x) 

Se f{r i )f(x j )< 0, então r l+1 =r, e s i+1 =x, 

Se f(r) i f(x j )> 0, então r M =x ; e s i+1 =s i 










Convergência 

Podemos observar que no Método da Bisseção determinamos uma raiz 
da equação, construindo sequências de intervalos ç e s, e uma sequência de 
soluções aproximadas x„/ = 0,1,... 

Essa sequência de soluções aproximadas é convergente para a solução 
desejada, uma vez que os intervalos são divididos pelos pontos médios corres¬ 
pondentes e são renomeados de forma que a raiz permaneça dentro do interva¬ 
lo. Ver mais detalhes na Referência [1]. 

Estimativa do número de iterações 

O número de iterações necessárias para se obter uma raiz x da equação 
f(x) = 0, pelo Método da Bisseção com uma precisão e > 0, previamente fixada, 
decorre do seguinte: 

Supondo-se que x está entre x n e s n , temos: 


X n- X \í( X n 



(s 0 -r 0 ) 

£n+1 


Impondo que 



para garantir que \x n -x 


<e temos: 


log 



\ 

< log (£) 


ou ainda, 


n> log(s 0 -r 0 )-log(e) ? 
log 2 

Logo, n é o número mínimo de iterações que devem ser realizadas para 
obter x com uma precisão e. 

Algoritmo 3.1 

1. Dados £>0, o intervalo inicial [ r 0 ,s 0 ] que contenha a raiz, isto é, 
f{r 0 )f(s 0 ) < 0. Faça Pare = Falso, / = 0 

2. Enquanto Pare = Falso faça: 

2.1 Determine x, = (ç +s),/2 











2.2 Se | f(Xj) | < £, então Pare = Verdade 
Senão 

Se f{r j )f{x i )< 0, então r i+1 =r) e s i+1 =x, 
Senão r M =x, e s l+1 = s, 

\x. — x • 

2.3 Se J—^—r-^- < £ < e então Pare = Verdade 

M 

Senão i = i + 1 


Exemplo 3.3 

Considere a seguinte função f(x) = 4x- e x . 

Usando o Método da Bisseção, podemos determinar x tal que f(x) = 0 
com uma precisão e = 0.01, conforme segue: 

a) Inicialmente determinamos graficamente uma vizinhança para a raiz, con¬ 
sideramos a forma equivalente 4x = e x , ou seja, f 1 (x) = 4x e f 2 (x) = e x , 
conforme ilustrado anteriormente na Figura 3.5. 

Observando a Figura 3.5, podemos concluir que a raiz x encontra-se na 
intersecção dos gráficos /j(x) e f 2 (x), e pertence ao intervalo [0, 1\ 

b) Considerando o intervalo inicial r 0 =0 e s 0 =1, temos f(0)f(1)<0, pois 
f(0) = -1 ef(1) = 1.2817 e, portanto, temos realmente uma raiz no inter¬ 
valo [0,1]. Observe que a função f(x) = 4x-e x é contínua no intervalo 
dado. 

c) Construindo a sequência de soluções aproximadas, temos que 
x 0 = (r 0 + s 0 ) /2 = (0 +1)/2 = 0.5000 é uma solução inicial. 

Como f(r 0 )f(x 0 ) = (-1)(0.3513)<0, temos que o novo intervalo ç e s 1 
será dado por: 

l r i= r o = ° 

}s, = x 0 = 0.5 


Calculamos a nova solução aproximada, conforme segue: 


x, = (ç +s 1 )/2 = (0 + 0.5 ) /2 = 0.2500 





Novamente, verificamos o critério de parada: 


x i~ x o 


= 1 > £, e como este 


não está satisfeito, repetimos o processo para calcular as novas soluções apro¬ 
ximadas, como segue: 

Como f(r 1 )f(x 1 ) = (-1)(-0.2840)>0 , temos que o novo intervalo r 2 e s 2 
será dado por: 


\r 2 =x 1 =0.25 
[s 2 = s r =0.5 


Calculamos o ponto médio do intervalo [r 2 , s 2 ] e temos: 
x 2 = (r 2 +s 2 )/2 = ( 0.25 + 0.5 ) /2 = 0.3750 


Verificamos novamente o critério de parada, dado por 



= 0.3333 > e, 


e como este não está satisfeito, sucessivamente repetimos o processo de cálculo 


das novas soluções aproximadas, como exibidas a seguir: 


x 3 = 0.3125 -» I X . 2 ' = 0.2000 > e 


x. = 0.3438 -» I * 3 ' = 0.0910 > e 

x 4 


x 5 = 0.3594 -> \ Xs * 4 I = 0.0434 > e 


x 6 = 0.3516 I Xe * 5 ' = 0.0222 > e 


x 7 = 0.3555 I X f Xg I =0.0110 > £ 


x T 


Xo = 0.3574 I X f X? I = 0.0053 < £ 

6 y 

*8 


Como o critério de parada está satisfeito, temos que a solução aproximada 
de f(X) = 0 é dada por x = x 8 = 0.3574. 
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Podemos calcular o número mínimo de iterações usando a expressão: 


log(s 0 -r 0 )-log( e) 
log 2 


1 = 5.6439 


Portanto, n > 5.6439, isto é, devemos executar no mínimo 6 iterações para 
obter a raiz x com a precisão c desejada, o que pode ser comprovado no exem¬ 
plo dado anteriormente com n = 7 iterações. 

3.3.2 Método de Newton 

Seja f: -»uma função contínua e diferenciável, desejamos determinar 
a raiz ou raízes de f(x), ou seja, determinar x tal que f(x) = 0. 

O Método de Newton consiste em, a partir de uma solução inicial aproxima¬ 
da x 0 dada inicialmente, determinar uma sequência {x,}, / = 1,2,... de soluções 
aproximadas para a raiz de f(x) da seguinte forma: 

Traçamos a reta tangente ao gráfico da função f(x) no ponto x,, e no pon¬ 
to em que essa tangente cortar o eixo das abscissas, temos a primeira solução 
aproximada x i+1 para a raiz x. 

Novamente no ponto x M traçamos a reta tangente ao gráfico de f(x), e no 
ponto em que essa nova tangente cortar o eixo das abscissas, temos a segunda 
solução aproximada para a raiz, e assim sucessivamente até determinarmos, 
com uma tolerância prefixada e, a raiz x desejada. 

Ilustramos graficamente, conforme Figura 3.9: 
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Figura 3.9 Método de Newton. 













Definindo como a o ângulo formado com o eixo das abscissas através da 
reta tangente ao gráfico da função f(x) no ponto x, (Figura 3.9), temos: 

\_ f(*i) 

J (*/-*/+») 

Do Cálculo Diferencial Integral, sabemos que tg(a) é a derivada da função 
f(x) no ponto x,, isto é, tg(a) = f\x i ) e, assim, podemos escrever: 



f\x i ) = 


f(Xj) 

(Xi~X i+1 ) 


Portanto, temos o processo iterativo chamado de Método de Newton, como 
segue: 


x i+ i=x i - 


f(Xj) 

f\x,) 


i = 0,1,2,... 


Observação 

O Método de Newton é também conhecido como Método das tangentes, 
devido a sua interpretação gráfica natural. 

Convergência do Método de Newton 

A convergência do Método de Newton pode ser tratada usando o Teorema 
de convergência do Método das Aproximações Sucessivas, o qual pode ser vis¬ 
to na Referência [1]. 

Definição 3.2 Convergência quadrática 

Dizemos que um método iterativo possui convergência quadrática se 

e i+ i 

lim—^- = k, em que k é chamada constante assintótica de proporcionali- 

e f 

dade, e j = | x ( - x e e j+1 = | x i+1 - x | são os erros cometidos nas iterações 
correspondentes. 

Teorema 3.1 

O Método de Newton possui convergência quadrática. 

Prova: Referência [1]. 

Para observar as boas propriedades de convergência do Método de 
Newton, exibimos o seguinte exemplo: 
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Exemplo 3.4 

O cálculo do número 43, usando o Método de Newton, consiste na resolu¬ 
ção da equação x 2 - 3 = 0 e, usando uma tolerância fixa com e = 10 4 , temos: 


*M = X i 


^1 = X, 

f (*;) 


(xf-3) 

(2x,) 


A partir de uma solução x 0 inicial, geramos a sequência de soluções 
aproximadas: 

x 0 = 1.5000 -» solução inicial dada 
x, = 1.7500 -> I X \ ~ X ° I = 0.1429 > e 

\ x i I 

x P = 1.7321 -» \ \ = 0.0103 >e 

| X 2 | 

X q = 1.7321 ^l = 0 . 0QQ0<e 

hfj 


Como o critério de parada está satisfeito, podemos parar e observar que 
essa sequência converge para x = 1.7321 = 43. 

Dessa forma, podemos observar que à medida que os valores de x k apro¬ 
ximam da raiz x, a convergência torna-se muito rápida, isso devido à proprieda¬ 
de da convergência quadrática do Método de Newton. 

Algoritmo 3.2 


1. Defina as funções f(x), f\x) e considere e > 0 uma tolerância fixa. 

2. Escolha x 0 uma solução inicial. 

Faça Pare = Falso e / = 0 

3. Enquanto Pare = Falso faça: 

f(Xi) 


3.1 x j+1 = x, - 


f (*,) 


\ X — X \ 

3.2 Se J— 4 —< e ’ então Pare = Verdade 

X M 
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Senão i = i + 1 











Observação 

O leitor pode incluir nesse algoritmo uma modificação no critério de para¬ 
da, considerando o valor da função no ponto x h isto é, | f(x ( .) | < e. 

Exemplo 3.5 

Usando o Método de Newton, resolvemos a equação cos(x) - x = 0 com 
£ = 0 . 001 . 


A partir do processo iterativo x i+1 = x, 
x 0 = 0.7 solução inicial dada 
x. = 0.7394 I X ! X ° I = 0.0533 > e 

\ x i I 

x ? = 0.7394 \ \ = 0.0000 > e 

\x 0 


f(x j ) „ . 

, ' , geramos a sequencia: 

f(x,) 


Como o critério de parada está satisfeito, tomamos como solução aproxi¬ 
mada para f{x) a solução x = x 2 = 0.7394. 

Observação 

Podemos, ainda, modificar o Método de Newton, da seguinte forma: 

O valor calculado da derivada na I a iteração, f\x 0 ) = k, em que ke % é 
fixado e substituído no processo iterativo de Newton durante as iterações: 

Assim, temos: 


o qual é conhecido como Método Modificado de Newton. 

3.3.3 Método das Secantes 


O Método das Secantes consiste em modificarmos o Método de Newton, 


cujo processo iterativo é dado por x j+1 = x, - 
função f'(Xj) da seguinte forma: 


f(Xj) 

Hx,) 


aproximando a derivada da 









f\x,) = 


(Xl-Xt-l) 


Observe que nesse caso, estamos trocando a inclinação da reta tangente 
pela inclinação da reta secante à curva, conforme pode ser observado Figura 
3.10 a seguir. 

Substituindo a expressão aproximada da derivada f\xj) no Método de 
Newton, podemos escrever: 

f(Xj) 

X — X — ___ !— _ 

1+1 1 fixj-fjx^) 

(*/”*/-») 


Simplificando a expressão anterior obtida, temos o Método das Secantes 
conforme segue: 

X ' +1 f(x i )-f(x j _ 1 ) 

Assim, dados os pontos x i _ 1 e x ( , traçamos a reta secante passando por 
(Xj_ v f{x i _ 1 )) e {x j ,f{x i )) \ onde esta cortar o eixo das abscissas temos a aproxi¬ 
mação x i+1 para a raiz x, conforme ilustrado na Figura 3.10. 



Figura 3.10 Método das Secantes. 

Assim, na Figura 3.10, considerando o ângulo a formado pela reta secante 
à curva no eixo das abscissas, podemos escrever: 
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ígK= '(*.-> = m-nx M ) 

(X,-X M ) I X - X, ,, 


Assim, temos o Método das Secantes dado por: 

_ Xj-i f(Xj) — x l f(Xj ? ) 

f(.Xj) — f(x i1 ) 

Algoritmo 3.3 

I.Seja f(x) contínua e e>0 uma tolerância fixa. 

2. Escolha x 0 e x 1 duas soluções aproximadas iniciais. 
Faça Pare = Falso e / = 0 

3. Enquanto Pare = Falso faça: 

31 x f(x l )-x j f(x l _ 1 ) 

1 + 1 f( X j ) — fi X j^ ) 

x. — x ■ \ 

3.2 Se J—< £, então Pare = Verdade 

I */+7 I 

Senão / = / + / 


Convergência 

Como o Método das Secantes é uma modificação do Método de Newton, 
as condições de convergência são parecidas, observando-se que não temos 
mais a propriedade de convergência quadrática. 

Quando f(x k ) = , podemos ter problemas de convergência, isto é, a 

sequência gerada pelo método pode divergir. 

Exemplo 3.6 

Considere a equação x 2 + 4x + 3 = 0] usando o Método das Secantes com 
e = 0.01, podemos resolver essa equação conforme segue: 

Podemos escrever a equação dada na forma equivalente x 2 =-4x-3. 

Chamando f 1 (x) = x 2 e f 2 (x) = -4x -3, temos uma vizinhança para as 
raízes, na intersecção dos gráficos de /j(x) e f 2 {><), conforme Figura 3.11: 










Figura 3.11 Localização de raízes. 


Usando o processo iterativo do Método das Secantes 
y x jf( x i~i) e tomando como soluções iniciais x n =-0.5 e 

x 1 =-2.0, obtemos a solução aproximada para a raiz, como segue: 


Xgfjx^-x^jXg) _ (-0.5)(-1)-(-2)(1.25) _ 3 

f( Xl )-f(x 0 ) (-1) - {1.25) -2.25 


Calculamos o critério de parada, dado por: 



0.5000 > £ 
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Como o critério de parada não está satisfeito, determinamos as demais 
soluções aproximadas, como seguem: 


x 3 = -0.5001 -» I X , X . 2 ' = 1.6661 > £ 


x, = -1.0769 -» I = 0.5356>e 


x 5 = -1.0159 -» I X ? X . 4 ' = 0.0600 > £ 


x 6 = -0.9996 -» I ^ Xg I = 0.0163 > e 


Xy = -1.0000 -> I X | Xg I = 0.0004 < £ 
















Como o critério de parada está satisfeito, temos a solução da equação 
x = -1.0000. 

Observação 

Podemos também determinar a outra raiz negativa, conforme exibida no grá¬ 
fico da Figura 3.10, apenas mudando as soluções iniciais e novamente usando o 
processo iterativo do Método das Secantes, até o critério de parada ser satisfeito. 

3.4 Exercícios 

1. Usando os métodos da Bisseção, de Newton e das Secantes, resolva a 
equação 3cos(x)-e °- 4x =0 e e = 0.1. Localize graficamente uma vizinhança 
para as raízes. 

2. Seja a função f(x) = e x ~ 2 +x 5 - 4. Determine o valor de x em que a fun¬ 
ção f(x) = 2, usando o Método de Newton com e = 0.001. 

3. Determine x tal que f(x) = 0 , em que f(x) = x 3 - cos(x) + 2, usando: 

a) Método da Bisseção com e = 0.1 

b) Método de Newton com e= 0.001 

c) Método das Secantes com e= 0.01. 

4. Determine um ponto da função f(x) = 6x 4 - 3x 3 -2x-3, em que a pri¬ 
meira derivada de f(x) se anula, usando o Método de Newton com e =0.001. 
O ponto encontrado é um ponto de mínimo da função dada? Justifique teorica¬ 
mente suas afirmações. 

5. Considere a função f(x) = 3x 2 + sen(x) - 2. Localize graficamente uma 
vizinhança para as raízes de f(x). 

a) Determine a raiz negativa de f(x) usando o Método de Newton com 
£ = 0 . 001 . 

b) Determine a raiz positiva usando o Método das Secantes, com 
£ = 0 . 1 . 

6. Determine as raízes da função f{x) = 0.2x 3 -3.006x 2 + 15.06x - 25.15, 
usando um dos métodos vistos com e = 0.001. 

7. Considere a função f(x) = x-sen(x)-0.5. Determine x tal que f(x) = 0, 
usando: 

a) Método das Secantes com e = 0.01. 




b) Método de Bisseção com e = 0.1. 

c) Determine uma raiz para f(x) dada, com o Software Numérico, usando 
o Método de Newton e e = 0.0001 - Referência [1]. 

8. A concentração de bactérias poluentes C em um lago diminui de acordo 
com a seguinte função: 


C(t) = 75e~ 15 ‘ + 20e-° ° 75t 


Deseja-se determinar o tempo necessário para que a concentração de 
bactérias seja reduzida ao valor de 15. 

a) Resolva o problema usando o método gráfico. 

b) Resolva o problema usando o Método de Newton com aproximação 
inicial t = 6e e = 0.001. 

9. Usando o Software Matlab, plote o gráfico e determine uma raiz das fun¬ 
ções abaixo: 

a) f(x) = 1 / [(x - 0.3) 2 + 0.01]+ [1/(x -0.9) 2 + 0.04]-6 

b) f(x) = e x -cos(x) -2 

c) f(x) = 2x 3 + ln(x) -5 

d) f(x) = x 5 -6.7x 4 + 8.4x 3 -10.8x 2 + 8x -6.8 

10. Faça um mapa conceituai detalhado sobre solução numérica de equa¬ 
ções, considerando os métodos vistos, introduzindo labeis indicando se aprendeu 
(A), não aprendeu (N), se gostou (G), detestou (D ), achou interessante (I), etc. 




UNIDADE 4 


Interpolação e Aproximação de funções 



4.1 Introdução 


Apresentamos nesta unidade a aproximação de uma função de uma va¬ 
riável real por outras funções mais simples, de modo que operações em geral 
sejam realizadas com mais facilidade. 

Essa situação ocorre quando temos uma função f(x) e esta apresenta 
um grau de dificuldade, por exemplo, para avaliar em pontos, derivar ou ainda 
integrar, ou mesmo quando conhecemos essa função em um número finito de 
pontos de um intervalo [a, b], geralmente obtida em experimentos, sem o co¬ 
nhecimento de sua forma analítica. 

Dessa forma, essa função será aproximada por funções polinomiais, expo¬ 
nenciais, trigonométricas, etc., que representarão a função original, e para obter 
qualquer informação sobre a função original, utilizamos a sua forma aproximada. 

Exibimos nesta unidade apenas a aproximação de uma função f(x), uti¬ 
lizando funções polinomiais, pela simplicidade no tratamento, continuidade e 
diferenciabilidade dessas funções, embora possamos estudar outros tipos de 
aproximações, como as indicadas nas Referências. 

Posteriormente, apresentaremos o Método dos Mínimos Quadrados para 
aproximação de uma função definida em um número finito de pontos, e nesse caso, 
consideramos a minimização da soma dos quadrados dos erros nos pontos. 

Como ilustração, apresentamos o problema da temperatura em um lago, 
conforme Figura 4.1. 



Figura 4.1 Lago. 

A temperatura T(°C) em um lago numa determinada região e época do ano 
é dada em função da profundidade P(m), conforme dados a seguir: 










T(°C) 

23.8 

22.8 

21.8 

20.8 

19.8 

18.7 

11.1 

P(m) 

0 

2.3 

4.9 

9.1 

13.7 

18.3 

22.9 


Deseja-se saber qual a profundidade do lago nas temperaturas 17.3°C e 
13.4°C, que não foram medidas. 

Como tais valores não se encontram na tabela, para resolver esse proble¬ 
ma, necessitamos do conhecimento de interpolação, como segue: 

4.2 Interpolação Polinomial 

Considere uma função f(x) definida em x 0 , x p x n , n + 1 pontos distintos 
de um intervalo [a, b ], e denotamos y, = f(x j ), i = 0,1,..., n conforme represen¬ 
tação na Figura 4.2: 



Figura 4.2 Valores da função f(x) nos pontos de interpolação. 

Desejamos aproximar a função f(x) definida em x 0 , x p ..., x n , n + 1 pontos 
distintos de um intervalo [a, b] por um polinómio P(x) de grau menor ou igual 
n, tal que este coincida com a função nesses pontos, isto é, P(x / ) = f(x / ) = y / , 
i = 0,1,...,n. 

O Teorema seguinte garante a existência e a unicidade do polinómio que 
desejamos determinar. 

Teorema 4.1 Existência e Unicidade 

Considere uma função f(x) definida em x 0 , x p ..., x n , n +1 pontos distintos 
de um intervalo [a, b], então existe um único polinómio P(x) de grau menor ou 
igual a n, tal que P(x,) = f(x ; ) = y„ i = 1,...,n. 





















Prova: 


Seja um polinómio de grau n, P(x) = a n x n + a n _ 1 x n 1 +.... + a 1 x + a 0 , tal 
que P(x,) = f(x j ) = y„ i = 0,1,..., n. 

Dessa forma, para que o polinómio coincida com a função nos n + 1 pon¬ 
tos, temos: 

' a n X 0 + a n-1 X V + 

a n X 1 + a n-1 X V + 

< 

ax n n + a„ iX 1 }- 1 + ... + a.x + a n = y n 

l n n n-1 n 1 n 0 J n 


+ a,x 0 + a 0 = y 0 
+ a 1 x 1 +a 0 = y 1 


Podemos observar que obtemos um sistema de equações lineares Ax = b, 
em que x = (a n , a n _„..., aj, b = (y 0 ,y 1 ,...,yj e a matriz A dada por: 


*3 

v n-1 
* 0 

.. X, 

x ? 

x ; n -' . 

.. X. 

X n 

v n-1 

.. X, 


O det(A), chamado de determinante de Vandermonde, é dado por 

det{A) = Y[{x j -Xj) = (x 1 -x 2 ){x 1 -x 3 )...{x n _ 1 -x n ) 

i<i 


Como os pontos x,-,/ = 0,1,...,n são distintos, podemos mostrar que 
det(A)#0, o que significa que o sistema linear possui uma única solução e, 
portanto, os coeficientes a 0 ,a p ...,a n do polinómio são únicos, dados pela reso¬ 
lução desse sistema. Assim, o polinómio P(x) existe e é único. 

Definição 4.1 Polinómio interpolador 

Denominamos polinómio interpolador de uma função f(x) definida em 
x 0 ,x p ...,x n , n+1 pontos distintos de um intervalo [a, b] o polinómio P(x) de 
grau menor ou igual a n, tal que coincide com a função nos pontos x (J / = n 
isto é, P(x j ) = f(x j ) = y j , i = 1,...,n. 

Representamos graficamente, conforme Figura 4.3: 








Figura 4.3 Interpolação. 

Embora o polinómio interpolador P(x) coincida com a função nos pontos 
de interpolação x 0 , x p ..., x„, espera-se que P(x) = f(x) para x*x,, i = 0,...,n, 
ou seja, estimamos f(x) pelo polinómio interpolador e cometemos um erro nessa 
aproximação, dado por: 

E(x) = f(x)-P(x) 

Podemos representar, graficamente conforme Figura 4.4: 



Figura 4.4 Erro de interpolação. 

Apresentamos a seguir uma expressão geral para o erro cometido quando 
aproximamos uma função f(x) por um polinómio interpolador P(x). 






















Teorema 4.2 


Seja f(x) uma função definida em x 0 ,x p ...,x n , n + 1 pontos distintos de 
um intervalo [a, b] e n + 1 vezes diferenciável. Se P(x) interpola f(x) nesses 
pontos, então o erro cometido em um ponto x é dado por: 

E(x) = f(x)-P(x) = -^~r + 1 \E 1 ) 

(n + 1)! 

n 

em que \\i(x) = Y\(x-x i ) e [x 0 , x n ]. 

i=0 

Prova: Referência [1]. 

4.3 Limitante superior para o erro 

Na expressão do erro do Teorema 4.2, o parâmetro £, não é conhecido no in¬ 
tervalo [x 0 , x n ] e, portanto, não é possível calcular o valor numérico de f {n+1) (^). 

Dessa forma, apresentamos uma estimativa para o erro como segue: 

E(x) = f(x)-P(x) = -^-f^ n+1 \E > ) 

(n + 1)! 

Podemos escrever: 


E(x) 


¥(*) 

(n + 1)! 


f (n+1) (E,) 


| V(x) | 

(n + 1)! 


f {n+1 \t,)\< l ¥(X) l M 
1 (n + 1)! 


com M= máx{| T n+? >(x) |, xe [x 0 ,x n ]|. 


Assim, temos um limitante superior para o erro: 


E(x) I < I y W i M. 


(n + 1)! 












Observações 

Podemos calcular uma estimativa para o erro somente quando tivermos a 
expressão analítica da função f(x), pois de acordo com a fórmula do limitante 
superior para o erro, temos que dispor da (n +1) -ésima derivada dessa função. 

Nos casos em que tivermos apenas a função tabelada em um número finito 
de pontos, sabemos que estamos cometendo um erro no ponto a ser avaliado, 
mas não é possível estimá-lo. 

Embora a resolução do sistema linear obtido na prova do Teorema 4.1 for¬ 
neça uma maneira para determinar o polinómio interpolador de uma função, apre¬ 
sentamos, também nesta unidade, a Fórmula Interpolatória de Newton-Gregory. 
Nas Referências podemos ver alguns estudos sobre outras fórmulas interpolató- 
rias, como a Fórmula de Lagrange e a Fórmula de Newton, as quais são fórmu¬ 
las diferentes, mas que levam no mesmo polinómio interpolador, pois como foi 
demonstrado, este é único. 

Para o desenvolvimento do polinómio interpolador de Newton-Gregory, 
apresentamos o conceito de diferenças finitas, como segue: 

4.4 Diferenças finitas 


Considere umafunção f(x) contínua no intervalo [a, b]. Sejam x 0 , x p ..., x n , 
n +1 pontos distintos deste intervalo [a, b ], tais que x j+1 - x, = h, i = 0,1,... n-1, 
isto é, os pontos são equidistantes. 

Definição 4.2 Diferença finita de ordem zero 

A diferença finita de ordem zero de uma função f(x) definida nos pontos 
x e [a, b] é dada por: 

A°f(x) = f(x) 

Definição 4.3 Diferença finita de ordem n 

A diferença dividida de ordem n de uma função f(x) definida nos pontos 
xe [a, b] é dada por: 


A n f(x) = A n 1 f(x + h) - A n 1 f(x) 

Nessa definição, para r = 1 temos o operador diferença finita progressivo 
A dado por: 
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A f(x) = f(x + h)-f(x) 


Assim, desenvolvendo-se os operadores podemos escrever 

A °f(x) = f(x) 

A 1 f(x) = f(x + h)-f(x ) 

A 2 f(x) = f(x + 2h) - 2f(x + h) + f(x) 

A 3 f(x) = f(x + 3h) - 3f(x + 2h ) + 3f(x + h)- f(x) 







í n l 


A n f(x) = 

f) 

f(x + nh)~ 


f(x + (n-1)h) + ... + (-1) n 


f(x) 


) 





Exemplo 4.1 

Considere uma função tabelada nos pontos, como segue: 


x < 

0.1 

0.2 

0.3 

f(Xi ) 

6 

12 

15 


Cálculo das diferenças finitas de ordem zero: 

A °f(x 0 ) = f(x 0 ) = 6 
A °f(x 1 ) = f(x 1 ) = 12 
A °f{x 2 ) = f{x 2 ) = 15 

Cálculo das diferenças finitas de ordem 1: 

a7(x 0 ) = A °f( Xl )- A °f(x 0 ) = (12-6) = 6 
A7(x ? ) = A °f(x 2 ) - A °f(x 1 ) = (15 -12) = 3 















Cálculo das diferenças finitas de ordem 2: 


A 2 f(x 0 ) = AV(x ; ) - a7(x 0 ) = (3 - 6) = - 3 

Podemos organizar o cálculo das diferenças finitas, conforme a Tabela 4.1 
a seguir: 

Tabela 4.1 Diferenças finitas. 

A°f A'f À 2 f À 3 f 


A°f(x 0 ) 

AT(x 0 ) 


A°f(x i ) 


A 2 f(x 0 ) 


ATfx^ 

A 3 f(x 0 ) 

A°f(x 2 ) 

AT(x 2 ) 

A 2 f (x i ) 

A°f(x 3 ) 




Os valores das diferenças finitas estão dispostos na tabela da seguinte forma: 

A 1 f(x 1 ) = A°f(x 1 ) - A°f(x 0 ) 

A 2 f(x 0 ) = À l f(x 1 )-À l f(x 0 ) 

A 3 f(x 0 ) = A 2 f(x 1 ) - Á 2 f(x 0 ) 

Exemplo 4.2 

Considere uma função f(x) tabelada nos pontos, como segue: 



0.2 

0.4 

0.6 

0.8 

f(x,) 

7.8 

5.9 

10.8 

13.2 


Conforme Tabela 4.1, construímos as diferenças finitas: 



M— 

CO 

<1 

H— 

OJ 

<1 

< 

O 

<1 

0.2 

7.8 


-1.90 

0.4 

5.9 6.80 


4.90 -9.30 

0.6 

10.8 -2.50 


2.40 

0.8 

13.2 
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4.5 Fórmulas Interpolatórias 

4.5.1 Fórmula Interpolatória de Newton-Gregory 


Considere uma função f(x) definida em um intervalo [a, b] e, x 0 , x p ..., x n , 
n + 1 pontos distintos e equidistantes deste intervalo. A Fórmula interpoladora 
de Newton-Gregory é dada pelo seguinte Teorema 4.3: 

Teorema 4.3 

Seja uma função f(x) definida em x 0 , x v ..., x n , n + 1 pontos distintos e 
equidistantes de um intervalo [a, b]. O polinómio 

_ 0 , A ’f(x n ) A 2 f(x n ) 

P(x) = A°f(x 0 ) + (x-x 0 ) + (x-x 0 )(x-x 1 ) 2jh 2 +- + 

+ (X - x 0 )(X - X,) ... (X - X n _-j) A ^ X h ° n ) 


interpola a função f{x) dada nestes pontos, isto é, P(x, ) = / = 0,1,... 

Prova: Referência [1]. 

Como desenvolvido anteriormente, uma estimativa para o erro é dada pelo 
limitante superior: 




| V(x) | 

(n + 1)! 


max 


f (n+1 \x) 


,xe [x 0 , x n ] 


Exemplo 4.3 

Considere a função f(x) = ln(x) tabelada nos pontos, como segue: 


x i 

2 

3 

4 

f( x i ) 

0.6931 

1.0986 

1.3863 


Podemos determinar o polinómio interpolador Newton-Gregory, avaliar 
f(3.2) e calcular um limitante superior para o erro, como segue: 

Usamos os três pontos tabelados para obter um polinómio interpolador de 
grau <2: 


P(x) = A°f(x 0 ) + (x-x 0 ) 


A V(x 0 ) 
Uh 1 


+ (x 


x 0 )( x 


*i) 


A 2 f(x 0 ) 

2 !h 2 
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Conforme Tabela 4.1, construímos as diferenças finitas: 



H— 

OJ 

<1 

< 

H— 

O 

<1 

2 

0.6931 


0.4055 

3 

1.0986 -0.1178 


0.2877 

4 

1.3863 


Assim, temos o polinómio interpolador de Newton-Gregory: 


P(x) = 0.6931 + (x-2) 0-4055 + (x - 2)(x - 3) 0/1178 


=-0.0589x 2 + 0.7000X -0.4713 


Portanto, f(3.2) = P(3.2)=1.1656 
Limitante Superior para o erro: 

I vi/íx) I , 

£(x) <-- l -máx A n+? >(x) xe[x n ,x 1 

1 (n + 1)! 1 1 0 " 


Como a função f(x) = ln(x), temos que as suas derivadas: 

fO)(x) = —,x*0, f (2 \x) = -—, x * 0 e A 3 >(x) = —. 
x x 2 x 3 


2 

Podemos observar que a função f (3 \x) = — é decrescente em módulo no 

x J 

intervalo [2,4], portanto, máx\f 3 (x)\ = 0.2500, em x = 2. em x = 2. 

Assim, 


E(3.2)< 


(■ 3.2-2){3.2-3){3.2-4) 
3! 


(0.25) = 0.0080. 


4.5.2 Polinómio Interpolador de Newton-Gregory com mudança de variável 


Podemos, ainda, escrever o polinómio interpolador de Newton-Gregory, 
considerando a seguinte mudança de variável: 












u = ——— ou x = x n +uh 
h 0 

Como os pontos são equidistantes, x r =x 0 + rh segue que: 

(x-x r ) = (u-r)h 

Note que nessa nova variável u os pontos de interpolação x 0 ,x 1 ,...,x n 
correspondem aos pontos 

u = 0,1,2,..., n. 

Dessa forma, podemos escrever o polinómio de Newton-Gregory na variável u: 

A 1 f(x n ) A 2 f(x n ) 

P{u) = A °f(x 0 ) + (u-0 ) — + (u- 0)(u -1) 2! ° +... 

A n f(x n ) 

+(U - 0){u — Y 1 

n! 


O limitante superior para o erro é dado por: 


E(u)< 


(n + 1)! 


h n+1 M 


em que M = máx\ f (n+1 \x) , xe [a,b]. 

Exemplo 4.4 

Considere a função f(x) da tabela nos pontos, conforme segue: 


*/ 

0.1 

0.2 

0.3 

f(Xi) 

2 

8 

10 


Podemos determinar o polinómio interpolador de Newton-Gregory na variá¬ 
vel u e avaliar f(0.12) como segue: 

O ponto x = 0.12 corresponde na variável u: 


u = 


(X-Xp) 


( 0 . 12 - 0 . 1 ) C2 


h 


0.1 
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. AV(x 0 ) _ A 2 f(x n ) 

P(u) = A°f(x 0 ) + (u -0) —^ + {u - 0)(u -1) —AyAZ 

Conforme Tabela 4.1, construímos a tabela de diferenças finitas: 



A°f A'f A 2 f 

0.1 

2 


6 

0.2 

8 -4 


2 

0.3 

10 


Assim temos: 


P{u) = 2 + u{6) + u{u-1)-^- = -2u 2 +8u + 2 


Portanto, f{0.12) = P u {0.2) = 3.52. 


4.6 Aproximação de funções: Método dos Mínimos Quadrados 

O Método dos Mínimos Quadrados consiste em, dados os valores de uma 
função f(x) nos pontos x„ i = m , determinar uma função g(x) que melhor 
se aproxime dessa função conhecida. Para essa aproximação, usaremos fun¬ 
ções polinomiais, exponenciais, logarítmicas, trigonométricas, etc. 

Assim, conhecidos os valores de f(x) nos pontos x jt i = 1,..., m, desejamos 
determinar g(x), tal que g(x j ) = f(x i ), i = 

Como ilustração, apresentamos o problema de previsão de escoamento de 
água, como segue, conforme Figura 4.5: 









Figura 4.5 Escoamento de água. 

A estimativa dos tamanhos dos reservatórios depende de estimativas acu¬ 
radas do escoamento da água no rio que está sendo confinado. 

Para alguns rios, registros históricos de longa duração de tais dados de 
escoamento são difíceis de obter. Entretanto, os dados meteorológicos sobre 
precipitação geralmente estão disponíveis mesmo após muitos anos. 

Portanto, é útil determinar a relação existente entre o escoamento e a pre¬ 
cipitação. Essa relação pode ser usada para fazer uma estimativa do escoamen¬ 
to nos anos em que apenas as medidas de precipitação foram feitas. 

Os dados disponíveis para um rio que deve ser estancado estão dispostos 
na seguinte tabela: 


Precipitação (cm) 

88.9 

108.5 

104.1 

139.7 

127 

94 

116.8 

99.1 

Escoamento (m 3 /s) 

14.6 

16.7 

15.3 

23.2 

19.5 

16.1 

18.1 

16.6 


Deseja-se determinar uma previsão do escoamento anual de água se a 
precipitação for de 120 cm. 

Caso discreto 

Consideremos o caso em que a função f(x) é definida em um conjunto 
discreto, isto é, a função é conhecida em m pontos, geralmente obtidos em expe¬ 
rimentos, conforme a tabela: 


x i 


x 2 



f(x,) 

f(*f) 

f(x 2 ) 


f(x m ) 





















Graficamente, temos a disposição dos pontos obtidos no experimento, con¬ 
forme Figura 4.6: 



Figura 4.6 Dados de um experimento. 

Observando a disposição dos pontos (x /; f(x,)), / = 1,..., m na Figura 4.6, 
vemos que g(x) possui o comportamento de uma reta, isto é, um polinómio de 
grau 1: 


g(x) = a 7 g 7 (x) + a 2 g 2 (x) = a,x + a 2 


Com g 7 (x) = x e g 2 (x) = 1. Assim, escolhemos uma família de funções as 
quais dependem dos parâmetros a 1 e a 2 . 

O problema agora consiste em determinar os parâmetros a 1 e a 2 de modo 
que a função g{x) melhor se ajuste aos dados da tabela. 

Para falar em melhor ajuste, temos que ter um critério para a escolha dos 
parâmetros a 1 e a 2 , isto é, temos que ter uma medida para o erro cometido 
nessa aproximação. 

Definição 4.4 Erro ou desvio 

Definimos e(x, ) = f(x,) — g(x, ) como o erro ou desvio, cometido numa apro¬ 
ximação de uma função f{x) por uma função g(x), nos pontos x ; , / = 1,..., m. 

Dessa forma, desejamos determinar uma função g(x) de modo que nos pon¬ 
tos x„ i = 1,..., m os desvios sejam pequenos. Nesse caso é tentador desejar que 

m 

a soma dos erros sejam mínimos, isto é, que X e ( x /) se J a m í nima - Entretanto, 

i=i 


102 









esse fato não traduz que g(x) seja uma boa aproximação para a função f(x), pois 

m 

podemos ter (x,) = 0 , sem que os erros sejam pequenos. 

i=i 

m 

Poderíamos, também, considerar ]T|e(x ( .)| m í nima ’ porém esse critério 

i=i 

acarreta dificuldades de resolução, pois a função valor absoluto não é diferenciá- 
vel na origem. 

Uma maneira para contornar esses problemas consiste em considerar uma 
medida para o erro da seguinte forma: 

m m 

Minimizar ^ e(x,) 2 = Minimizar (f(x,) -ç/(x,)) 2 

i=i i=i 

Assim, considerando o exemplo da Figura 4.6, desejamos encontrar uma 
função g(x) = a 1 x + a 2 que melhor se aproxime da função f{x), de forma que 

m 

E(a p a 2 ) = ^(e(x, )) 2 seja mínimo. Do Cálculo Diferencial de funções de várias 

i=i 

variáveis, se a função E{a p a 2 ) possui um ponto de mínimo, então suas deriva¬ 
das parciais devem ser nulas nesse ponto, isto é: 


dE 

da 1 


= 0 


dE n 

e - = 0 

9a, 


Derivando E(a v a 2 ) com relação à variável a 1 , temos: 


9 E_ = _d_ 
9a, 9a, 


i=i 


9a, 


£(e(x,.)) 2 =— £(a,x, +a 2 -f(x i )) í 


í=i 


= X 2{a 1 x j + a 2 -f(x j ))x j 


i=i 


= 2 


a,£xf+a 2 £x ; 


i=1 


i=1 


-2 


5>M> 


i=1 


= 0 


Assim, temos: 


f m 2 ^ 


( m 

X*/ 

a, + 




l /=í ) 


i=l 


i) 

















Derivando E(a p a 2 ) com relação à variável a 2 , temos: 


ÓE 

da 2 


c) 

da 2 



(e(x ; .)) 2 


\ 


v 


d 

da 2 


r m 




+ ^2 


y 


= X 2 (¥í + 


i =1 


a 2 


f(Xi)) 


= 2 


m m 

a /S x / +m(a 2 )-^f(x i ) 


i =1 


i =1 


= 0 


Assim, temos; 

m 

a 1 +ma 2 = Y J f{x i ) 

í=i 

Portanto, os parâmetros a, e a 2 que minimizam o erro E(a p a 2 ) necessa¬ 
riamente satisfazem o seguinte sistema de equações lineares: 



< 


f m 2 N 


' m ^ 

X*/ 

a, + 

x*, 

J 



r m ^ 



l x i 

a 1 + ma 2 = 


a 2 = X 


X f (*/) 

i =1 


Xi f( x i ) 


O sistema de equações obtido é chamado sistema de equações normais, o 
qual pode ser resolvido por qualquer método numérico visto anteriormente. 

Exemplo 4.5 

Considere a seguinte função f(x) da tabela como segue: 


x i 

0.5 

0.8 

0.9 

1.0 

1.1 

f(x,) 

1.51 

1.82 

1.91 

2.00 

2.10 


Plote os pontos dados e verifique se os mesmos possuem um comporta¬ 
mento linear, isto é, de uma reta (verifique isso!). 

Usando o método dos Mínimos Quadrados, determinamos dentre todas as 
retas g(x) = a r x + a 2 aquela que melhor se ajusta aos dados. 





























5 5 

I*, 2 X*, 

i=1 i=1 


1 


X /(x /K 

i=1 

1? 

\x 

Oi 

1_ 


A. 


-1 

_1 


Temos a seguinte tabela de dados: 



*1 

f(Xi ) 

x? 

f(x i )x j 

z 

0.5 

1.51 

0.25 

0.7550 

0.8 

1.82 

0.64 

1.4560 

0.9 

1.91 

0.81 

1.7190 

1.0 

2.00 

1.00 

2.0000 

1.1 

2.10 

1.21 

2.3100 

4.30 

9.34 

3.91 

8.2400 


Assim, temos o seguinte sistema de equações lineares: 


'3.91 

4.30 ' 

V 


8.24 

4.30 

5 

. a 2. 


9.34_ 


Usando o Método de Eliminação de Gauss, na forma pivotada, temos: 

3.91 4.30 . 8.24 

0 0.2711 . 0.2785_ 

Assim, 

a 2 =1.0272 
a 1 = 0.9778 


Portanto, a função linear que melhor se aproxima dos dados da tabela é 
dada por g(x) = 0.9778x + 1.0272. 

Cálculo do erro: 

e(x,) 2 = (f(0)-g(0)) 2 =0.0000 
e{x 2 ) 2 = {f(1)-g(1)) 2 =0.0001 
e(x 3 ) 2 = (f(2) - g(2)) 2 = 0.0000 
e(x 4 ) 2 = (f (3)-g(3)) 2 =0.0001 
e(x 5 ) 2 = (f (4) -g(4)) 2 =0.0000 
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5 

Portanto, e(x,) 2 = 0.0001 e qualquer outra reta possui a soma dos qua- 

i =1 

drados dos erros superior a esse valor obtido. 

Podemos ter dados experimentais em que seja necessário aproximar a 
função f(x) por um polinómio de grau 2, isto é, uma parábola: 

g(x) = a l g 1 (x) + a 2 g 2 {x) + a 3 g 3 (x ) = a r x 2 +a 2 x + a 3 


com g 1 (x) = x 2 ,g 2 (x) = x eg 3 (x) = 1. 

Generalizando esse procedimento, escrevemos g(x) como uma combina¬ 
ção linear de funções como segue: 


g(x) = a^fx) + a 2 g 2 {x) +... + a n g n {x) 


com as funções g,(x) escolhidas. 

Procedendo de maneira análoga ao caso do ajuste linear, a reta, podemos de- 

m 

terminar os parâmetros a j ,i = 1,...,n de forma que o erro ^e(x,) 2 seja mínimo. 

i=i 

Assim, temos: 


dE < m 
— = 0 o 
da 1 


i=i 


2 \9Mi)9Mi) a ; + J j g 2 (x i )g 1 (x j )) 


+ 


y v /=í 

A m 


ã 2 + . 


i =1 




i=l 


dE _ ( m 

■ = 0 <=> 


da 0 


i =1 


J, 9 1 (x i )g 2 (x i ) a 1 + ^g 2 (x,)g 2 (x,) 


+ 


y v /=í 

A m 


ã 2 + . 


i=1 


a n='L f ( x i )g 2 ( x i) 


í=i 


dE _ ' m 

- = 0 o 


da„ 


\ í 


)$„(*,) a i + 


i=i 


i =1 


y v' 

f m \ m 

5> n (*,&„(*,) a n =X^/)a„(>f/) 


a 2 +... 


í=i 


y 


i=i 


+ 
























Portanto, para determinar os parâmetros a„i = 1 ...,n temos que resolver 
o seguinte sistema de equações lineares: 


A C 


'Zg 1 ( x i)g 1 ( x i) a,+ J\g 2 (x i )g 1 (x i ) 


y v 


a 2 +... 


f m 


'Lg n (x l )g 1 ( x i) a n = 'Lf(x i )g 1 (x i ) 

i=1 i=1 


f rn 


'Zg 1 ( x i)g 2 ( x i) 


f m 


y Zg 2 ( x i)g 2 ( x i) 


Z gJx, )g 2 ( x i ) a n = £ f ( x i)g 2 ( x i ) 


i=1 


i=1 


t,gi( x i)g n ( x i) 

v'= 7 y 


a i + 


í m \ 

Tg 2 ( x i)g n ( x i') 


a 2 +... 


í m 


a n = Z^ x /)5- n ( x /) 

i=1 i=1 


O sistema de equações lineares obtido é denominado sistema de equações 
normais, o qual pode ser resolvido por qualquer método visto anteriormente. 

Resolvido esse sistema, determinamos os parâmetros a,, / = e, con¬ 
sequentemente, a função: 


g(x) = a 1 g 1 (x) + ... + a n g n (x) 


que melhor se ajusta à função f{x) nos pontos x p ...,x m no sentido dos míni¬ 
mos quadrados. 

Exemplo 4.6 

Considere uma função f(x) definida conforme tabela: 


*/ 

-2 

-1 

0 

1 

2 

3 

f( x i) 

19.01 

3.99 

-1.00 

4.01 

18.99 

45.00 


Dispondo os pontos num gráfico podemos ver que g{x) pode ser ajustada 
por uma parábola. Verifique esta afirmação! 

Assim, tomamos g(x) = a 1 x 2 +a 2 x + a 3 , isto é, g 1 (x) = x 2 , g 2 (x) = x, 
g 3 (x) = 1 e determinamos os parâmetros a p a 2 e a 3 de modo que g(x) se ajus¬ 
te aos dados da tabela no senso dos mínimos quadrados. 
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Temos o sistema de equações normais: 


6 6 6 

xv 1 .x, 3 x*, 2 

1=1 i=1 i=1 

fí fí fí 


V 


ÍHx,)xf 

i=1 

Q 

O O O 

É */ 3 É x / 2 É*/ 

i=1 i=1 i=1 


a 2 

= 

È f ( x ,) x / 

i=1 

hl 2 í*i 6 

_ i=i i=i 


a 3_ 


S^) 

. i=i 


Temos a seguinte tabela de dados: 




x? 

x? 

xf 

f(Xi) 

x, f(Xi) 

xff(Xi) 


-2 

4 

-8 

16 

19.01 

-38.02 

76.04 


-1 

1 

-1 

1 

3.99 

-3.99 

3.99 


0 

0 

0 

0 

-1.00 

0.00 

0.00 


1 

1 

1 

1 

4.01 

4.01 

4.01 


2 

4 

8 

16 

18.99 

37.98 

75.96 


3 

9 

27 

81 

45.00 

135.00 

405.00 

3 

19 

27 

115 

90.00 

134.98 

565.00 


Assim, podemos escrever o sistema de equações normais, como segue: 


~115 

27 

19 

V 


565.00 

27 

19 

3 

ã 2 

= 

134.98 

19 

3 

6 

a 3_ 


90.00 


Usando o Método de Eliminação de Gauss, na forma pivotada, temos: 

~115 27 19 . 565.00 

0 12.6604 -1.4612 . 134.98 

0 0 2.6926 . 90.00 

Assim: 

a 1 = 5.0893 
< a 2 = 0.0515 
a 3 =-1.1403 
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Portanto, g(x) = 5.0893x 2 + 0.0515 x-1.1403 


































Cálculo do erro: 


e(x r ) 2 = (f(x 1 ) —g(x 1 )) 2 =0.0108 
e(x 2 ) 2 = (f(x 2 ) -g(x 2 )) 2 = 0.0086 
e(x 3 ) 2 =(f(x 3 ) -g(x 3 )) 2 =0.0197 
e(x 4 ) 2 = (f(x 4 ) -g(x 4 )) 2 = 0.0001 
e(x 5 ) 2 = (f(x 5 ) -g(x 5 )) 2 = 0.1088 
e(x 6 ) 2 = ( f(x 6 ) -g(x 6 )) 2 = 0.0332 

6 

Portanto, e(x, ) 2 =0.1812 e qualquer outra parábola ajustada possuem 

i=i 

a soma dos quadrados dos erros superiores a esse valor obtido. 

4.7 Exercícios 

1. Seja uma função f(x) definida conforme a tabela abaixo: 


X 

0 

1 

2 

3 

f(x) 

-5.6 

-4 

8.1 

9 


Determine o polinómio interpolador de f(x), avalie f(2.3) usando a Fórmula 
Interpolatória de Newton-Gregory. 

2. Seja uma função f(x) definida conforme a tabela abaixo: 


X 

0.1 

0.2 

0.3 

0.4 

f(x) 

8 

3.9 

4.8 

-8.9 


Determine o polinómio interpolador de f(x) nos pontos dados, usando a 
Fórmula Interpolatória de Newton-Gregory na variável u. Avalie f(0.35). 

3. Estude o polinómio interpolador de Lagrange e o limitante superior para 
o erro - Referência [1]. 

4. Seja a função f(x) = e x + cos(x) +1 tabelada como segue: 


X 

0 

0.5 

1.1 

f(x) 

3 

3.53 

4.46 


Usando o polinómio interpolador de Lagrange, avalie f(0,6) e um limitante 


superior para o erro. 
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5.0 calor específico da água em função da temperatura em °C é dado por: 


t 

30 

35 

40 

45 

c 

0.99826 

0.99818 

0.99828 

0.99890 


Com o Software Numérico - Referência [1] usando o polinómio interpola- 
dor de Newton-Gregory, calcule aproximadamente o calor específico para 

t = 37.5 °C. 

6. A partir dos seus conhecimentos de aproximação de funções, resolva o 
problema da temperatura em relação à profundidade de um lago, proposto no 
início deste capítulo, usando um polinómio interpolador de grau 2. Analise os 
resultados obtidos. 

7. Usando o Método dos Mínimos Quadrados, determine g(x) = a 1 x + a 2 
que melhor se ajusta aos dados da tabela abaixo: 


X 

-2 

-1 

0 

1 

2 

3 

4 

f(x) 

-4.98 

-3.00 

-1.01 

0.99 

3.01 

4.98 

7.01 


Calcule ^e(x ( ) 2 e analise os resultados obtidos. 

i=1 

8. Com o Software Numérico - Referência [1 ] —, e usando o Método dos Mí¬ 
nimos Quadrados, determine uma função g(x) que melhor se ajusta aos dados 
da tabela abaixo: 


X 

-2 

-1 

0 

0.1 

1 

2 

3 

f(x) 

6.00 

1.01 

0 

0.12 

3.00 

9.98 

20.98 


7 

Calcule ^,e(x ( ) 2 e analise os resultados obtidos. 

i=i 

9. Considere o problema de previsão de escoamento de água, exibido neste 
capítulo: 

a) Faça um gráfico dos dados da tabela. 

b) Usando o Método dos Mínimos Quadrados, determine a melhor reta 
que melhor se ajusta aos dados da tabela e faça uma previsão do es¬ 
coamento anual de água se a precipitação for de 120 cm. 

10. Usando o Software Matlab, resolva os exercícios 7) e 8). 




























11 .Faça um Mapa Conceituai detalhado sobre os tópicos de interpolação, 
introduzindo labeis indicando se aprendeu (A), não aprendeu (N), se gostou (G), 
detestou (D), achou interessante (I), etc. 

12. Faça um Mapa Conceituai detalhado sobre os tópicos de ajuste de cur¬ 
vas - MMQ introduzindo labeis indicando se aprendeu (A), não aprendeu (N), 
se gostou (G), detestou (D), achou interessante (I), etc. 





UNIDADE 5 


Integração Numérica 



5.1 Introdução 


Nesta unidade, apresentamos alguns métodos numéricos para calcular 
aproximadamente o valor da integral de uma função com uma variável real defi¬ 
nida num intervalo fechado [a, b]. 

De maneira geral, temos: 

b 

1 = J f(x)dx 

a 

Em que a função f(x) é contínua com derivadas contínuas no intervalo 
[a, b\ 

Sabemos que o valor da integral I é dado pelo Teorema Fundamental do 
Cálculo para integrais como: 

b 

I = J f(x)dx = F(b) - F(a) 

a 

em que F(x), chamada de função primitiva de f(x) é tal que F\x) = f(x). 

b 

Graficamente, podemos representar I = J f(x)dx como a área A, formada 

a 

entre a curva e o eixo das abscissas quando f(x) assume valores positivos em 
[a, b], conforme a Figura 5.1 a seguir: 



Figura 5.1 Área sob o gráfico da função f(x). 

Exemplo 5.1 

' x 6 , 

Para calcular \x dx, como F(x) = — satisfaz F(x) = x , então: 
i 6 
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J x 5 dx = F(1) - F(0 ) = —- = 1/6 

o 6 6 

Muitas vezes, o cálculo da função primitiva torna-se algo trabalhoso, em al¬ 
guns casos conhecemos apenas os valores da função tabelados em um número 
finito de pontos e não podemos calcular a integral dessa função usando apenas 
os conhecimentos do Cálculo Diferencial Integral, como citamos anteriormente. 

Métodos numéricos são desenvolvidos para calcular aproximadamente o 

b 

valor da integral I = J f{x)dx, em que é necessário apenas o conhecimento da 

a 

função f{x) em um número finito de pontos. Esse processo é conhecido como 

Fórmulas de Quadratura de Newton-Cotes. 

Para ilustrar o tópico de integração numérica, consideremos o problema de 
estimar medidas das margens até um rio conforme Figura 5.2: 



A partir de uma linha reta próxima às margens de um rio, foram feitas me¬ 
didas em (m) entre essa linha reta e as margens do rio, de 15 em 15 metros, a 
partir do ponto tomado como origem. Tais dados foram registrados conforme a 
tabela abaixo. 



0 

15 

30 

45 

60 

M, 

0 

1 

0 

2 

4 

m 2 

5 

7 

7 

8 

10 


Deseja-se determinar o valor aproximado da área coberta pelo rio. 















5.2 Integração numérica usando interpolação 


Este método consiste em aproximar uma função f{x) definida em 
x 0 ,x p ...,x n , n + 1 pontos distintos do intervalo [a, b], pelo polinómio interpola- 
dor P n (x), e integrar o polinómio P n (x) ao invés da função f(x). 

Assim, temos: 

X X 

J f (x)dx = J P n (x)dx 


Podemos representar graficamente, conforme Figura 5.3: 



Figura 5.3 Integração numérica. 


5.3 Fórmulas de Quadratura de Newton-Cotes 


Considere uma função definida em x 0 , x p ..., x n , n + 1 pontos distintos e 
equidistantes do intervalo [a,b], isto é, x j+1 - x, = h, i = 0,...,n-1, com h>0 a 
distância entre os pontos. 

Nesse caso, podemos considerar o polinómio interpoladorde Newton-Gregory 
para f(x)\ 


P(x)=D°f(x 0 ) + (x-x 0 ) 


DVjx,,) 

1! h 


+ (x-x 0 )(x-x 1 ) 


D 2 f(x 0 ) 
2! h 2 


+... + (x — x Q )... (x — x n1 ) 


A n f(x n ) 


n!h n 













Dessa forma, temos: 


jf(x)dx= J P n (x)dx 


Como f(x)=P n (x) + E n (x) em que o erro é dado por: 


( n+1) 

EJx) = (x-x„)...(x-x n )riAL-, x 0 <%<x 


Temos 


J f{x)dx = J P n (x)dx + J E n (x)dx 


Como os pontos são distintos e equidistantes, podemos considerar o poli¬ 
nómio de Newton-Gregory na variável u, em que: 

u = (x - x 0 )/h 
Assim, 


P„{u) = A°f(x 0 ) + uA 1 f(x 0 ) + U(U 2! — A 2 f(x 0 ) +. 

+ u(u-1)...(u-(n-1)) 

n! °' 


e o erro na interpolação na variável u é dado por: 


= u(u-1)...(u-n) 

" (n+1)! 0 " 


Dessa forma, o erro na integração pode ser escrito como: 


E n = \ E n {x)dx = h] E n (u)du = h\ 


u{u-1)...{u-ri) 
(n + 1)! 


(n+1) 


h n+1 f(x) du 









Erro cometido na Integração Numérica 


X, 


X, 


Podemos aproximar J f(x)dx por j P n (x)dx e nessa aproximação come- 

x o x o 

temos um erro expresso pelo seguinte teorema: 

Teorema 5.1 

a) Se f(x) possui n +1 derivadas contínuas no intervalo [ x 0 , x n ], e os pon¬ 
tos x, = x 0 + ih, i = 0,1,...,n subdividem o intervalo de integração num 
número ímpar de intervalos, então existe um ponto t ,, tal que: 



0 


b)Se f(x) possui n + 2 derivadas contínuas no intervalo [x 0 , xj, e os pon¬ 
tos x, =x 0 + ih, i = 0,1,...,n subdividem o intervalo de integração num 
número par de intervalos, então existe um ponto E ,, tal que: 



Prova: Referência [1]. 

5.3.1 Regra dos Trapézios 

Considere uma função f(x) definida em dois pontos x 0 e x 1 no intervalo 
[a, b]. O polinómio interpolador da função f(x) degrau n = 1 é dado por: 



e 


X, 


J f{x)dx = J p 1 (x)dx = hjp 1 (u)du 


X, 


0 


em que u = ——^ 
H h 
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Representamos graficamente, conforme Figura 5.4: 



Figura 5.4 Regra dos Trapézios. 


Assim, 


"7 ' 

J f(x)dx = hj (A °f{x 0 ) + ívAY(x 0 )) du 


= hj A °f{x 0 )du + h j i/A 7 f(x 0 )du 


Como A°f(x 0 ) = f{x 0 ) e A 1 f{x 0 ) = f{x 0 )-f(x 0 ), temos 


7 ' ' 

J f(x)dx = /ij A°f{x 0 )du + hj uA 1 f(x 0 )du 


= hf(x 0 )u 


h(/(x ; )-/(x 0 )) 


u 




: hf(Xg ) + — (f(Xj ) — f(Xg )) 


:-(f(x 0 ) + f(x,)) 



























Portanto, temos: 


j f(x)dx = -T [f(x 0 ) + f{x 1 )] 


denominada Regra dos Trapézios. 

Observação 

O cálculo da integral J f{X)dx corresponde à área do trapézio formada 

x 0 

entre o polinómio interpolador P ? (x) e o eixo das abscissas, conforme a Figura 
5.4, justificando, assim, a denominação Regra dos Trapézios. 

Erro na Regra dos Trapézios 

Nesse caso, como o intervalo de integração n = 1 é ímpar, e usamos a parte 
(a) do Teorema 5.1, temos: 


E i = ^jí—j(u)(u-1)du, x 0 <í<x 1 


i 

Como ju(u-1)du 

0 


—, temos que o erro é dado por: 


E,= 


12 


f (2) & 


Limitante superior para o erro 

Podemos observar que o argumento £, na fórmula do erro, não é uma 
grandeza numérica conhecida no intervalo [x 0 , x y ] e, portanto, não é possível 
calcular o erro cometido na Regra dos Trapézios. Dessa forma, podemos calcu¬ 
lar uma estimativa para o erro, que será o limitante superior para o erro: 


Ei 


12 




12 


[X 0 ,x ; ], 


Como | P 2) (^) | < máx{| f {2) (x)\, xe [x 0 ,x r ]j-, temos que: 


E 1 |<^máx{| P 2 >(x) |,x 0 <x<x ? } 
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Exemplo 5.2 


2 

Podemos calcular J (sen(x) + x)dx usando a Regra dos Trapézios e um limi- 

i 

tante superior para o erro da seguinte forma: 

Tabelando a função f(x) nos pontos x 0 =1.0 e x 1 = 2, temos: 


X 

1 

2 

f(x,) 

1.8415 

2.9093 


Usando a Regra dos Trapézios, temos: 

2 h 1 

J (sen(x) + x)dx = — [f(x 0 ) + f(x ? )] = —[1.8415 + 2.9093] = 2.3754 


Limitante superior para o erro: 

E 1 \ < — máx{\ A 2, (x)|> 7 < x < 2} 


( 2 ) 

Como a função f(x) = -sen(x), então 


f(2) 



= máx {| T 2 >(x) | ,1<x<2} = 1 


pois a função | f (2 \x) | = sen(x) assume o valor máximo em x = 

Dessa forma, temos: 

i i (D 3 

E. <^-(/) = 0.0833 

5.3.2 Regra dos Trapézios Generalizada 

A Regra dos Trapézios Generalizada consiste em dividirmos o intervalo 
de integração [a, b] em n subintervalos iguais, de amplitude h, de forma que 
x 0 =a e x n =b e, aplicarmos a Regra dos Trapézios em cada um deles, isto é, 
em cada dois pontos consecutivos. 
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Graficamente, temos: 



Figura 5.5 Regra dos Trapézios Generalizada. 

x " h h 

J f(x)dx = -[f(x 0 ) + f(x 1 )] + - \f(x 1 ) + f(x 2 )] + ... 


+^lf(x n - 1 ) + f(x n )l = 

= — [f(x 0 ) + 2f(x 1 ) + 2f(x s ) + ... + 2f(x n _ 1 ) + f(x n )] 


Erro total na Regra dos Trapézios Generalizada 

O erro total na fórmula generalizada é obtido a partir da soma dos erros 
cometidos em cada subintervalo, isto é: 

Como vimos em cada aplicação da Regra dos Trapézios, temos a seguinte 
expressão para o erro: 

E i=-j^ f ^x 0 <^< Xl 


Assim, o erro total cometido é dado por: 

n tj3 (2) 

/=í ' 

Como x j _ 1 < ç ( . <Xj, i = 1, n e f {2 \x ) é uma função contínua por hipóte¬ 
se, então existe £ e [x 0 ,x n ], tal que 


123 









X f (2) (^/) = nf (2) (£)* X 0 -Ç>-X n 

i=1 

Assim, a expressão para o erro na Regra dos Trapézios Generalizada 
torna-se: 

h 3 (2) 

E T=-jj nf (&^[ x o’ x nl 

Como o número de subintervalos n é dado por n = (x n -x 0 )/h, temos: 

h 2 (2) 

E t =-—(x n - x 0 )f(^),^e[x 0 , x n ] 

Limitante superior para o erro 

Como o argumento £ não é conhecido, não podemos determinar o erro 
exatamente, mas podemos calcular um limitante superior para o erro. 

t?2 

| E t \ ^J^( X n- X o) máx {\ fi 2 ) (X)\,X 0 <X<X n } 


Exemplo 5.3 

Usando a Regra dos Trapézios Generalizada, podemos calcular o valor 

i 

aproximado da integral J(x/n(7 + x))dx com 5 subintervalos e um limitante su¬ 
ei 

perior para o erro. 

Solução: 


n = 



—> 5 = 


( 1 - 0 ) 

h 


->h = 0.2 


Dessa forma, podemos tabelar a função f(x) como segue: 


*1 

0 

0.2 

0.4 

0.6 

0.8 

1.0 

f(Xi) 

0 

0.0365 

0.1346 

0.2820 

0.4702 

0.6931 














Assim, temos: 


7 f_ 

J (x ln(1 + x))dx= — [f(x 0 ) + 2(f(x 1 ) + f(x 2 ) + f(x 3 ) + f(x 4 )) + f(x 5 )] = 

= ^-(2.5397) = 0.2538 

Limitante Superior para o erro: 

| E t \^( X n- X o) máX {\ f(2) ( X )\’ X 0^ X ^ X n} 

(2 + x} 

Como a função f {2) = --/ é uma função decrescente em módulo em 

( 1+x) 2 

[0 1] (verifique isso, observando que a função f {3 \x) < 0 no intervalo de integra¬ 
ção), então máx\ A 2 Hx) | = 2 em x = 0. 

Assim, 

\E t \< ^(x n - x 0 ) máx{| A 2 >(x) |, 0 < x < l}= °£(1 - 0)(2) = 0.0333 

5.3.3 Regra 1/3 de Simpson 

Considere uma função f(x) definida em três pontos distintos x 0 , x 1 e x 2 
e equidistantes no intervalo [a, b] e o polinómio interpolador da função f(x), de 
grau n = 2, com x 0 = a e x 2 = b. 

Assim, 

_ A 7(x„) A 2 f(x n ) 

P 2 (x) = A°f(x 0 ) + (x-x 0 ) 0 +(x-x 0 )(x-x r ) ° 

e 

*2 X2 2 

J f(x)dx = j P 2 (x)dx = hjP 2 {u)du 

x o x o 0 


em que u = 


(x-x 0 ) 


h 
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Graficamente, temos: 



Figura 5.6 Regra 1/3 de Simpson. 
Assim, 


J f{x)dx = hj í (A °f(x 0 ) + uA 1 f(x 0 ) + U(U f 7) A 2 f{x 0 ) 

x 0 O V * * 

2 2 

= hj A °f{x 0 )du + hj uA 1 f(x 0 )du 

O 0 

+ h 2 j^_H & 2 f(Xo)du 

O * ■ 

Sabendo-se que 

Af(x 0 ) = f(x 0 ), A J f(x 0 ) = f(x 1 ) - f(x Q ) 

A 2 f(x 0 ) = f(x 2 ) + 2f(x 1 ) + f(x 0 ) 

Temos: 

x 2 2 2 2 

J f(x)dx = h j A° f(x 0 )du + h J uA V(x 0 )dtv +/? J 


du = 




= hf{x 0 ){u) 


A?(/(x,)-/(x 0 )) 






0 






















L. 

- [f(x 0 ) + 4f(x 1 ) + f(x 2 )] 

O 


Assim, obtemos: 


J f(x)dx = ^[f(x 0 ) + 4f(x 1 ) + f(x 2 )] 


denominada Regra 1/3 de Simpson. 

Erro na Regra 1/3 de Simpson 

Nesse caso, o intervalo de integração foi subdividido em um número par 
de subintervalos, n = 2k, k = 1, 2,3,... e portanto, pela parte (b) do Teorema 5.1, 
temos: 

( n+2) 



0 


Assim, para n = 2 temos: 



o 



0 



Como o argumento Ç não é conhecido, não é possível calcular o erro exa¬ 
tamente, portanto, trabalhamos com um limitante superior para o erro: 
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Limitante Superior para o erro 

Temos que o erro é dado por: 


e 2 = - 


( 4 ) 

h 5 f& 

90 


<x<x„ 


Tomando em módulo a expressão dada anteriormente, temos que: 


E 2 = 


h*_ 

90 


f i4) &) Ue [x 0 ,x 2 \ 


Dessa forma, podemos escrever o limitante superior para o erro como: 


h 5 


E 2 |<—{máx| n 4 \x) \,x 0 <x<x 2 } 


Exemplo 5.4 

2 

Usando a Regra 1/3 de Simpson, podemos calcular je x dx e um limitante 
superior para o erro da seguinte forma: 0 

Tabelando a função f(x) com h = 1.0, temos: 


x i 

0 

1.0 

2.0 

f( x i) 

1 

2.7183 

7.3891 


Dessa forma, temos: 

2 

J f(x) dx = - [f(x 0 ) + 4f{x 1 ) + f(x 2 )] 

0 J 

= -[1 + 4(2.7183 ) + 7.3891] = 6.4208 
3 

Portanto, 


2 

J e x dx = 6.4208 

o 












Limitante superior para o erro 

Da fórmula do limitante superior para erro, temos que: 



h5 

<— \máx I A 4 >(x) I 
90 1 1 1 


xe[0,2 ]} 


Como a função f (4) (x) = e x , então 


| f {4) (2.0) | = máx{| A 4) (x)|, x e [0,2]}= 7.3897, 

pois a função | f (4) (x) | é crescente no intervalo [0,2], 
Assim, 

I I <^7.3891 = 0.0821 
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5.3.4 Regra 1/3 de Simpson Generalizada 

A Regra 1/3 de Simpson Generalizada consiste em dividirmos o intervalo 
[a, b ] de integração em n subintervalos iguais de amplitude h, em que n é um 
número par de subintervalos, de forma que x 0 = a, x n = b , e aplicarmos a Regra 
1/3 de Simpson em cada 2 subintervalos, isto é, a cada 3 pontos consecutivos. 

Graficamente, temos: 



Assim, aplicando a Regra 1/3 de Simpson para cada 3 pontos, isto é, a 
cada 2 subintervalos, temos: 















j f(x)dx = ~ [f(x 0 ) + 4f(x 1 ) + f(x 2 )] + ^ \f(x 2 ) + 4f(x 3 ) + f(x 4 )] + ... 

*o J J 

+^[f(x n . 2 ) + 4f(x n _i) + f(x n )] 

= | [f(x 0 ) + 4(f(x 1 ) + f{x 3 ) + ... + f(x n1 )) 

+2{f{x 2 ) + f(x 4 ) + ...+ f(x n _ 2 ) + f(x n )] 

Erro total na Regra 1/3 de Simpson Generalizada 

O erro total cometido na Regra 1/3 de Simpson Generalizada é obtido a 
partir da soma dos erros cometidos a cada aplicação da Regra 1/3 de Simpson 
em dois subintervalos subsequentes, isto é: 

Para cada dois subintervalos, 3 pontos distintos e equidistantes, aplicamos 
a Regra 1/3 de Simpson e o erro cometido para cada aplicação é dado por: 

( 4) 

h 5 f(L) 

E 2 = - x 2i _ 2 < Ç, < x 2i , i = 1,2,..., n/2 

Assim, temos: 

n/2 h 5 (4) 

f( Ái )> x 2,-2 * ti Z x 2i, i = 1,2,...,n/2 

í=i yu 

Como x 0 <'%! < x n , i = 1,2,..., n/2 e a função f {4 \x ) é contínua por hipóte¬ 
se, então existe um [x 0 , x n \ , tal que: 

n/2 (4) 

)= n fm,x 0 <^<x n 


n/2 (4) 


Assim, a expressão para o Erro na Regra 1/3 de Simpson Generalizada é 
dada por: 


h 5 n <4) 

Et= ~9Õ^ m ’ X o-^- Xn 






(x — X ) 

Como o número de subintervalos n = — n 0 , temos que: 


M ( 4) 

E f =-—(x n -x 0 )f(fy,x 0 <í<x n 
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Limitante Superior para o erro 

Como [x 0 ,x n ] não é uma grandeza numérica conhecida, trabalhamos 
com um limitante para o erro, conforme segue: 

h 4 

I Et l“7sõ (Xn_x ° )máx {l f(4)(x) \ ,x o- x - x "} 


Exemplo 5.5 

Usando a Regra 1/3 de Simpson Generalizada, 

f 1 

J- dx com 4 subintervalos e um limitante superior para 

Q (l + X) 

forma: 


podemos 
o erro, da 


calcular 

seguinte 


Temos h = (Xf1 Xp) =- = 0.5 
n 4 


Tabelando a função f(x) com h = 0.5, temos: 


Xi 

0 

0.5 

1.0 

1.5 

2.0 

f(Xi) 

1 

0.6666 

0.5000 

0.4000 

0.3333 


Temos 

2 j #_ 

J - <lx= - [f(x 0 ) + 4f(x 1 ) + 2f(x 2 ) + 4f{x 3 ) + f(x 4 )] = 

0 (' + X) d 

= —[1 + 4((0.6666 ) + 0.4000) 

3 

+2(0.5000) + 0.33333] = 1.1000 


Assim, 


J 


_/_ 

(1 + x) 


dx = 1.1000 


















Limitante Superior para o erro 

Temos 


E t\ -Tfíf) {Xn ~ Xo) máx i\ fi4)(X) I’ x ° - x - x "} 
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24 

Como r 4) (x) = -- é uma função decrescente em módulo no intervalo 

(7 + x) 5 

[0,2], máx\ f< 4 \x) | = 24, em x=0. 

Assim, 

0 5 4 

\E t < ^-{2 - 0)(24) = 0.0166 
I M 180 

5.4 Exercícios 


0.7 

1. Calcule I = J ( e- 6x + 8.4x + 3)dx usando: 

0.1 

a) Regra dos Trapézios com 6 pontos. 

b) Regra 1/3 de Simpson com 6 subintervalos. 

c) Calcule um limitante superior para o erro em cada caso. 


2. Seja uma função f(x) tabelada como segue: 


X 

-1 

0 

1 

2 

3 

4 

5 

f(x) 

1 

0.51 

0.42 

0.82 

1.91 

0.99 

1.88 


5 

Calcule J f(x)dx usando a Regra dos Trapézios e a Regra 1/3 de Simpson. 

i 

3. Calcule j(4.6x 2 +cos(x)) + 3.5dx usando a Regra dos Trapézios com 6 

o 

pontos e um limitante superior para o erro. 

4. Estude a Regra 3/8 de Simpson e um limitante superior para o erro e 
resolva o exercício 2) - Referência [1]. 

2 

5. Calcular J( e 3x +8.5x + 6)dx usando: 

i 


a) Regra dos Trapézios com 4 subintervalos. 

















b) Regra 1/3 de Simpson com 6 subintervalos. 

c) Regra 3/8 de Simpson com 9 subintervalos e um limitante superior para 
o erro. 

6. Determine o menor número de subintervalos em que podemos dividir o 
intervalo [0 1] para obter 

i 

I = j (5 e- 2x )dx 
0 

usando a Regra dos Trapézios com erro menor ou igual a 0.001. 

7. Resolva o exercício 6) usando a Regra 1/3 de Simpson, a Regra 3/8 de 
Simpson e para a divisão encontrada calcular a integral I. 


8. Uma linha reta foi traçada de modo a tangenciar as margens de um rio 
nos pontos A e B. Para medir a área entre o rio e a reta AB, foram traçadas per¬ 
pendiculares em relação a AB com um intervalo de 0.05 m. Qual é essa área? 


Perpendiculares 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Comprimento(m) 

3.28 

4.02 

4.64 

5.26 

4.98 

3.62 

3.82 

3.68 

3.26 


Faça uma interpretação gráfica do problema e resolva-o usando seus co¬ 
nhecimentos de Cálculo Numérico. 

9. Com o Software Matlab, usando a Regra dos Trapézios, calcule as inte¬ 
grais das seguintes funções: 

5 

a) j(sen(x)-x 2 + 3.1)dx 

i 

3 

b) j (e x +cos(x)- 2.8)dx 

10. Com o Software Numérico - Referência [1] -, usando a Regra 1/3 de 
Simpson, calcule o valor aproximado da integral das seguintes funções: 

4 

a) J( e x + 5x + 4.9)dx , usando 6 subintervalos e um limitante superior 

i 

para o erro. 

2.5 

b) J (x 3 +6x + 2)dx , usando 9 pontos e um limitante superior para o 


erro. 















11.Faça um Mapa Conceituai detalhado sobre métodos numéricos para 
Integração Numérica, introduzindo labeis indicando se aprendeu (A), não apren¬ 
deu (N), se gostou (G), detestou (D ), achou interessante (I), etc. 




UNIDADE 6 


Solução Numérica de Equações Diferenciais 
Ordinárias e Parciais 



6.1 Introdução 


O estudo de Equações Diferenciais deve-se ao fato de que, em situações 
práticas nas diversas áreas científicas, frequentemente temos problemas cuja 
modelagem matemática é representada por uma equação diferencial. 

O desenvolvimento de métodos numéricos para resolução dessas equa¬ 
ções diferenciais é necessário, pois soluções analíticas nem sempre são possí¬ 
veis de serem determinadas. 

Tipicamente, problemas que envolvem uma função de uma variável e suas 
derivadas levam a uma Equação Diferencial Ordinária, e problemas que envol¬ 
vem uma função de várias variáveis e suas derivadas parciais nos levam a uma 
equação diferencial parcial. 

Um exemplo simples de equação diferencial ordinária é o problema da pro¬ 
pagação da fumaça de cigarro num ambiente fechado ou semifechado. 

Exemplo 6.1 

Considere uma sala de trabalho que contenha um volume V = 36 m 3 de ar 
originalmente isento de monóxido de carbono. Suponha que, a partir do instante 
t = 0, a fumaça do cigarro contendo 4% de monóxido de carbono é introduzida 
na sala por um fumante, com uma vazão de 0.002 m 3 por minuto, e que essa 
sala possui uma janela ou porta semiaberta, de forma que a mistura gasosa, 
distribuída homogeneamente, sai da sala com a mesma vazão. 

a) Qual é a expressão para a quantidade q(t) de monóxido de carbono na 
sala num instante t? 

b) Sabendo-se que a exposição prolongada de monóxido de carbono com 
concentração acima de 0.00012 é prejudicial à saúde, em quanto tempo 
essa concentração será atingida? 


Para equacionarmos o problema, consideremos a quantidade q(t) de mo¬ 
nóxido de carbono na sala no instante t. A taxa de variação de q(t) é igual a quan¬ 
tidade de monóxido de carbono que está entrando menos a que está saindo. 


dq 

dt 


- 0.002 

100 


0.002 

36 


< 7(0 


A solução analítica desse problema com a condição inicial q(0) = 0 é a função 


200 200 - 

q(t) = -e 

139 139 


139 t 


2500000 
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Considerando = temos: 

36 


139 t 

f(t) = 0.04- 0.04e 2500000 

que é a porcentagem da quantidade de monóxido de carbono na sala no instante 
t, como exibido na Figura 6.1: 



Figura 6.1 Monóxido de carbono numa sala de fumantes. 

Para sabermos o instante em que a sala atinge o limite insuportável de 
monóxido de carbono, resolvemos a seguinte equação: 

f(t) = 5Êi = 0.00012 

36 

E obtemos o instante t = 54.1 minutos, o que significa que em 54 minutos 
a sala estará saturada de monóxido de carbono. 

Entretanto, muitas equações diferenciais não possuem solução analítica, o 
que nos leva ao desenvolvimento de métodos numéricos para obter a solução 
aproximada. 

Definição 6.1 

Uma Equação Diferencial Ordinária de ordem n é uma equação da seguinte 
forma: 


F(t,y(t),y'(t),y"(t),...,/ n \t)) =0 










em que estão envolvidas a função incógnita y = y(f) e suas derivadas até a 
ordem n, sendo t a variável independente. 

A notação y {>) representa a derivada de ordem j da função incógnita y em 

dW 

relação a variável independente t e pode também ser representada por 

Exemplo 6.1 

As equações a seguir podem ser escritas conforme a definição anterior. São 
equações diferenciais ordinárias, em que a função desconhecida é y = y(x)e 

dy d 2 v 

as suas derivadas são y' = —,y" =——. 

7 dx dx 2 

dy 

a) — = 2x - y é uma equação diferencial de ordem 1. 
dx 

b) y"+(1- y 2 )y '+ 6y = 0, em que y = y(t ) e suas derivadas y '(t) = — 

dt 

d 2 y 

e y"(t) = , e uma equação diferencial de ordem 2. 

6.2 Equação diferencial de primeira ordem com valor inicial 

Um problema de valor inicial consiste em uma equação diferencial 
y' = f(x,y), e uma condição inicial y (x 0 ) = y 0 , em que y 0 é um valor da função 
y dado no ponto x 0 . 

Assim, temos um problema de valor inicial denotado por: 

\y' = f(x,y) 

0 ) 

{y(x 0 )=y 0 

Resolver o problema (1) consiste em determinar uma função y(x) que 
satisfaça (1), isto é, y' = f(x,y) e y(x 0 ) = y 0 . 

Exemplo 6.2 

Considere o problema do crescimento da poluição numa cidade e o modelo 
matemático em que o crescimento da poluição está sujeito à hipótese de que a 
taxa de variação da poluição p(t) é proporcional à poluição num instante t, dado 
pela equação diferencial: 

^ = kp , ou p\t) = kp(t) 


em que k é a taxa de variação da poluição anual. 
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Suponha que no instante t = 0 temos uma poluição inicial de Ippm (uma 
parte por milhão). Se k = 0.065 por ano, podemos utilizar métodos analíticos ou 
numéricos adequados de resolução e fazer a previsão dessa poluição num perío¬ 
do de tempo t = 10 anos, resolvendo o problema de valor inicial: 


p' = 0.065 p 
P(0) = 1 


A solução analítica dessa equação diferencial é dada pela função 
p(f) = e 0,065t , pois ao substituirmos p(f) e sua derivada na equação dada, veri¬ 
ficamos que a equação está satisfeita. 

Uma condição que garante a existência e a unicidade da solução do proble¬ 
ma de valor inicial (1) é dada pelo seguinte Teorema 6.1: 

Teorema 6.1 

Considere uma função real f(x, y) contínua no intervalo a<x<b, com a e 
b finitos e — o° < y < + oo. Se existe uma constante L tal que para todo xe [a, b] 
e para todo par de valores (y, y ; ) tivermos 

| f(x,y)-f(x,y 1 ) |<Z_| y-y 1 \ 

conhecida como Condição de Lipschitz, então existe uma única função y = y(x) 
satisfazendo as seguintes condições: 

a) y(x) é contínua e diferenciável para todo xe [a,b ]; 

b) y' = f(x,y (x)) para xe[a,õ]; 

c ) y(* 0 ) = y 0 , em que y 0 é um valor conhecido. 

Prova: Referência [11] 

Observação 

Podemos observar pelo Teorema 6.1 que a solução de uma equação dife¬ 
rencial ordinária é uma função y = y(x) contínua e diferenciável que satisfaz a 
equação diferencial y’ = f(x,y) e passa pelo ponto (x 0 ,y 0 ). 




6.3 Discretizagão 


O processo de discretizagão consiste em resolver numericamente uma 
equação diferencial calculando aproximações para a função y = y(x) em pon¬ 
tos discretos x 0 ,x v x 2 ,...,x N do intervalo de discretizagão [a,b]. 

Para discretizar o intervalo [a, b], tomamos N subintervalos, N>1 , e faze- 

X — X 

mos x n =x 0 +nh, n = 0,1,2,..., N com x 0 =a e x N =b, sendo h= N N ° ■ 

Conhecendo y(x) num ponto inicial dado y(x 0 ) = y 0 (valor inicial), calcu¬ 
lamos passo a passo, nos pontos x 1 =x 0 + h, x 2 =x 0 +2h, x 3 =x 0 +3h,..., 
x n =x 0 +nh, as soluções aproximadas y n para a solução exata y(x n ), 
n = 0,1,2,...,N, conforme ilustrado na Figura 6.2: 



Figura 6.2 Discretizagão de um intervalo. 

O erro local, cometido nas aproximações em cada ponto, é a diferença en¬ 
tre o valor exato da equação diferencial e o valor numérico aproximado em cada 
um dos pontos do intervalo [a, b], isto é: 


e(x n ) = y(x n )-y n ,n = 1,...,N 



















6.4 Métodos Numéricos para Equações Diferenciais Ordinárias 

6.4.1 Método de Euler 

Seja a equação diferencial de primeira ordem com valor inicial: 
íy' = f(x,y) 

i y(x 0 )=y 0 

O Método de Euler é baseado na expansão da função derivável y(x) em 
Série de Taylor, isto é, expandindo a função y(x) por Série de Taylor nas vizi¬ 
nhanças do ponto x n , até a ordem 1, podemos escrever: 

y(x n + h) = y(x n ) + hy\ x n ) 

Dessa forma, como o ponto x n + h = x n+1 e y'(x n ) = f(x n ,y(x n )), temos: 
y(x n+1 ) = y(x n ) + hf(x n ,y(x n )) 


Assim, podemos escrever para n = 0: 


y(Xj) = y(x 0 ) + hf(x 0 , y(x 0 )) 

e 


y(Xi) = yo + hf(x 0 ,y 0 ) 


Chamando y 1 =y(x 0 ) + hf(x 0 , y(x 0 )), temos o valor aproximado da função 
y(x) no ponto x p isto é, y(x 1 ) = y 1 . 

Portanto, o Método de Euler no primeiro passo calcula y 1 =y 0 + hf(x 0 ,y 0 ) 
que é a aproximação y 1 da solução exata y(x r ) no ponto x f . 

Repetindo esse procedimento para o ponto x 2 , temos no segundo passo 
y 2 =y 1 + hf {x p y r ) que é a aproximação y 2 da solução exata y(x 2 ) no ponto 
x 2 - E assim sucessivamente para cada um dos pontos x n , n = 0,1,2,...,N de 
um intervalo [a,b]. 




De maneira geral, temos: 


yn + 1=yn +h ^ X n’Vn)’ n = 0, l 2,... 

Podemos interpretar graficamente o Método de Euler aproximando a fun¬ 
ção y(x) pela reta tangente ao gráfico de y no ponto dado ( x 0 , y(x 0 )), da se¬ 
guinte forma: 


(y-y 0 ) = m(x-x 0 ) 


em que m = y'(x 0 ) = f(x 0 ,y 0 ) é o coeficiente angular dessa reta. 
Assim, podemos escrever: 

( y-y 0 ) = f ( x o’yo)( x - x o ) y=y 0 + f(x 0 ,y 0 ){x-x 0 ) 


Tomando y = y 1 e x = x p temos que: 


Yi =yo +f ( x o’y 0 )( x i- x o) 


Ou ainda 

yi=y 0 +hf( x o,y 0 ) 

Assim, o ponto (x 1 ,y 1 ) pertence a essa reta tangente, e o mesmo raciocí¬ 
nio pode ser utilizado para os demais pontos discretizados, conforme ilustramos 
na Figura 6.3: 



Figura 6.3 Método de Euler. 
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O erro local cometido é a diferença entre o valor exato e o valor aproximado 
em cada ponto (. x n ,y n ),n = 0,1,...,N. 


Erro de Truncamento 


por: E = 


No Método de Euler, truncamos a Série de Taylor em p = 1, e o erro é dado 

y"m 2 


2 ! 


[a b\. 


Como o parâmetro Ç não é conhecido no intervalo de discretização, vamos 
considerar uma estimativa para o erro, isto é, um limitante superior para o erro: 


E< 


Mh 2 

2 ! 


M = máx {|y "(x)|, x e [a, d]) 


Algoritmo 6.1 

Considere a equação diferencial, com valor inicial 


í y' = f(x,y) 

\y(x 0 ) = y 0 


1. Declare 


a) Função f(x,y) 

b) Condições iniciais: y(x 0 ) = y 0 

c) Intervalo [a,b], em que a = x 0 

_ Q 

d) Número de subintervalos N e calcule h = - 

' N 

2. Para n = 0,...,N -1, faça: 
início 

x n+1 =x n + h 

y n+ i=yn + hf ( x n’y n ) 

fim 


Exemplo 6.3 

Considere a seguinte equação diferencial de ordem 1 com valor inicial: 


y' = f(x,y) = y 

y(x 0 ) = y(0) = 1 








a) Usando o Método de Euler, podemos calcular a solução aproximada 
y(1.0 ) com h = 0.2. 

b) Compare os valores exatos dessa equação, com os valores aproxima¬ 
dos obtidos em a). 

Considerando h = 0.2, temos a discretização do intervalo dada por: 
x 0 = 0, x 1 = 0.2, x 2 = 0.4, x 3 = 0.6, x 4 = 0.8 ex 5 =1.0 

Do Método de Euler: 

yn + 1 =yn +hf (Xn,y„) 

Cálculo de y, 

Para n = 0, temos a condição inicial y(x 0 ) = y(0) = 1, então: 

Yi = y 0 + hf (x 0 ,y 0 ) 

= y 0 + 0.2 (y 0 ) = 1 +0.2(1) = 1.2000 

Cálculo de y 2 

Para n = 1, temos: 

y 2 = yi + hf(x 1 , y ; ) 

= y, + 0.2 (y 1 ) = 1.2000 + 0.2(1.2000 ) = 1.4400 

Cálculo de y 3 

Para n = 2, temos: 

y 4 = y 3 +hf(x 3 ,y 3 ) 

= y 3 + 0.2(y 3 ) = 1.7208 + 0.2(1.7208) = 2.0736 

Cálculo de y 4 


Para n = 3, temos: 




y 4 = y 3 +hf(x 3 ,y 3 ) 


= y 3 + 0.2(y 3 ) = 1.7208 + 0.2 ( 1.7208 ) = 2.0736 


Cálculo de y 5 

Para n = 4, temos: 


y 5 = y 4 + hf(x 4 ,y 4 ) 

= y 4 + 0.2(y 4 ) = 2.0736 + 0.2(2.0736 ) = 2.4883 


Portanto, temos que y(1) = y 1 = 2.4883. 

Estimativa para o erro 


E<—máx {| y"(x) |,xe [0,/]} 


Como a solução analítica y(x) da equação é dada por y(x) = e x , temos que: 
máx{| y"(x) |,xe [0,1]\= 2.7183, em x = 1, uma vez que y"(x) = e x é uma 
função crescente em módulo. 

Portanto, temos: 

E < {0 ' 2 ^ (2.7183) = 0.0544. 


Exibimos na Tabela 6.1 os valores exatos, a partir da solução exata 
y(x) = e x , a solução numérica pelo Método de Euler e os respectivos erros lo¬ 
cais nos pontos discretos. 


Tabela 6.1 


N 

*/ 

Sol. Exata 

y(*n) 

Sol. Aprox. 

y n 

Erro 

\y(x n )-y n \ 

0 

0 

1.0000 

1.0000 

0.0000 

1 

0.2 

1.2214 

1.2000 

0.0214 

2 

0.4 

1.4918 

1.4400 

0.0518 

3 

0.6 

1.8221 

1.7280 

0.0941 

4 

0.8 

2.2255 

2.0736 

0.1519 

5 

1.0 

2.7183 

2.4883 

0.2300 
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Observe na Tabela 6.1 que temos uma boa estimativa para o erro no início 
do processo, mas piora quando os pontos afastam-se do ponto inicial. 

Podemos observar também que para valores de h menores que h = 0.2 , 
obteremos melhores resultados, pois o erro depende da amplitude do passo h. 


6.4.2 Método de Euler Aperfeiçoado 


O Método de Euler Aperfeiçoado é uma modificação do Método de Euler 
desenvolvido anteriormente, da seguinte forma: 

1 ^ 

y n+ i =y n + T,V( x n’y n ) +f ( x n +^ /„ + h y n h n= 0 , 1 , 2 ,... 

Para melhor entendimento da expressão do Método de Euler Aperfeiçoado, 
faremos uma interpretação gráfica, conforme Figura 6.4 a seguir: 



Figura 6.4 Método de Euler aperfeiçoado. 

Observando a Figura 6.4, vemos que: 

a) A reta t passa pelo ponto A = (x n ,y n ) e possui coeficiente angular 
y n = f(x n ,y n ). Usando o Método de Euler, calcula-se y^ +1 =y n +hy' n \ 

b) A reta r passa pelo ponto B = (x n+1 ,y^ +1 ) e possui coeficiente angular 

f(x n +h,y n +hy' n ) = f(x n+v y E n+1 y, 

c) A reta r 1 passa pelo ponto B e sua inclinação é a média das inclinações 
das retas ter; 


d) A reta r 2 passa pelo ponto A e é paralela a reta ry 
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e) O valor y n+1 obtido pelo Método de Euler Aperfeiçoado é uma aproxima¬ 
ção para a solução y(x) no ponto x n+1 . 


Algoritmo 6.2 

Considere a equação diferencial 


íy'(x) = f(x,y) 

ly(* 0 ) = y 0 


1. Declare: 


a) Função f{x, y) 

b) Condições iniciais: y(x 0 ) = y 0 

c) Intervalo [a, b], em que a = x 0 

d) Número de subintervalos N e calcule: h = - 

' N 

2.Para n = 0,1,2,3,...,N -1 

Calcule: 

início 

X n+ 1 = X n +h 
k 1 = f ( X n’V n ) 

k 2 =f( x n + i,y n + hk i) 

y n+ i =y n + ^ k i +k 2 ) 

fim 

Exemplo 6.4 

Considere a equação diferencial com valor inicial: 


y' = f(x,y) = y 


y(x 0 ) = y(0) = 1 


Podemos, usando o Método de Euler Aperfeiçoado, determinar a solução 
aproximada y{1.0) com h = 0.2 da seguinte forma: 

Considerando h = 0.2, temos a discretização do intervalo: 
x 0 = 0, x 1 = 0.2, x 2 = 0.4,x 3 = 0.6, x 4 = 0.8 ex 5 =1.0 






Cálculo de y, 

Para n = 0 temos a condição inicial y(x 0 ) = y(0) = 1, então 
k^fiXcYo) -+k 1 =y 0 = 1 

k 2 = f(x 0 + h,y 0 + hk 1 ) = y 0 + hk 1 =1 + 0.2(1) = 1.2000 

Do Método de Euler Aperfeiçoado temos: 

1 ^ 

yn + 1=yn+-^(k 1 + k 2 ) 

Portanto: 

y 1 = y 0 +^(k 1 +k 2 ) = 1 + ^~ (1 + 1.200) = 1.2200 

Cálculo de y 2 

Para n = 1 

k 1 =f(x 1 ,y 1 ) ->k 1 =y 1 =1.2200 

k 2 = f(x 1 + h,y 1 + hk 1 ) = (y 1 + hky) = 1.2200 + 0.2(1.2200) 

= 1.4640 

y 2 =y 1 +^(k 1 +k 2 ) = 1.2200+ ^-(1.2200+ 1.4640) = 1.4884 

Cálculo de y 3 

Para n = 2 

k 1 =f(x 2 ,y 2 ) -+k 1 =y 2 = 1.4884 

k 2 = f(x 2 +h,y 2 + hk 1 ) = (y 2 + hk 1 ) = 1.4884 + 0.2(1.4884) 

= 1.7861 

y 3 =y 2 +^(k 1 +k 2 ) = 1.4884 + (1.4884 + 1.7861) = 1.8159 

Cálculo de y 4 


Para n = 3 




k 1 = f ( x 3 ’Y 3 ) -» k 1 =y 3 = 1.8159 


k 2 =f( x 3 + h,y 3 + hk 1 ) = (y 3 + hk 1 ) = 1.8159 + 0.2(1.8159) 

= 2.1791 


y 4 = y 3 + t ^(k 1 +k 2 ) = 1.8159 + ^-(1.8159 + 2.1791) = 2.2154 

Cálculo de y 5 

Para n = 4 

k = f{x 4 ,y 4 ) -» k 1 =y 4 =2.2154 

k 2 =f(x 4 + h, y 4 +hk 1 ) = (y 4 + hk^ = 2.2154 + 0.2(2.2154) 

= 2.6585 

y 5 =y 4 + l l(k 1 +k 2 ) = 2.2154 + 2.2154 + 2.6585) 

= 2.7028 


Portanto, temos que y(1.0) = y 1 = 2.7028. 

Exibimos na Tabela 6.2 os valores exatos e os valores aproximados obtidos. 


Tabela 6.2 


N 

*i 

Sol. Exata 

Sol. Aprox. 

Erro 



y(*n) 

y n 

\y(* n )-y n \ 

0 

0 

1.0000 

1.0000 

0.0000 

1 

0.2 

1.2214 

1.2200 

0.0014 

2 

0.4 

1.4918 

1.4884 

0.0034 

3 

0.6 

1.8221 

1.8159 

0.0062 

4 

0.8 

2.2255 

2.2154 

0.0101 

5 

1.0 

2.7183 

2.7028 

0.0155 


Como podemos observar na Tabela 6.2, o Método de Euler Aperfeiçoado 
apresentou melhores resultados quando comparado com o Método de Euler, isto 
é, os erros são menores. 

Para obtermos melhores resultados devemos diminuir a amplitude do passo 
h, por exemplo, tome h = 0.01 e resolva o exercício da mesma forma, como an¬ 
teriormente, e compare os resultados obtidos. 
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Observação 

Os Métodos de Euler e Euler Aperfeiçoado são chamados de Métodos de 
Runge-Kuttade I a ordem e 2 a ordem, respectivamente. Esses métodos são bas¬ 
tante utilizados, pela sua simplicidade e pela precisão dos resultados. 

Existem também, os Métodos de Runge-Kutta de 3 a e 4 a ordem, os quais 
podem ser encontrados com detalhes na Referência [2], 


6.5 Solução Numérica de Equações Diferenciais Parciais: 
Uma Introdução 


Como no caso das equações diferenciais ordinárias, vários problemas nas 
diversas áreas científicas podem ser formulados por uma equação diferencial 
parcial. 


Em geral, as equações diferenciais parciais são mais complexas do que as 
equações diferenciais ordinárias, pois envolvem um número maior de variáveis. 

Assim, quando a função incógnita depende de mais de uma variável e 
relaciona com suas derivadas parciais através de uma equação, temos uma 
equação diferencial a derivadas parciais. 

Por exemplo, a equação — -2^-^ = 0 é uma equação diferencial parcial 

dt dx 

do tipo parabólica em que a função incógnita é u = u(x,t). 

Nesta seção, veremos um método numérico para resolução das equações 
diferenciais parciais, usando para isto uma discretização por diferenças finitas. 


Apresentaremos um método de diferenças finitas para resolução de uma 
equação diferencial parcial, o qual é baseado na discretização do domínio e a 
substituição das derivadas presentes na equação diferencial por aproximações 
envolvendo valores numéricos da função. 


Para isso, inicialmente tratamos a aproximação das derivadas de uma função 
de uma variável por diferenças finitas, usando a expansão da Série de Taylor. 

Considere uma função y(x) contínua com derivadas contínuas, podemos 
escrever a função y(x) nas vizinhanças de um ponto xpor: 


y(x + h) = y(x) + hy '(x) + 


2 


y"{><) + ■■■ 


Ou ainda, 


h 2 


y(x - h) = y(x) - hy\x) + —y"(x) + ... 
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Tomando a expansão da função y(x) por Série de Taylor, nas vizinhanças 

do ponto xaté a primeira ordem, e desprezando os termos de ordem 0(h 2 ) e 
de ordem maior, temos: 


y(x + h) = y(x) + hy '(x) 


Dessa forma temos a aproximação da derivada da função y(x), da seguinte 
forma: 


y'(x) 


y(x + h)-y(x) 
h 


Considerando Ay(x) = y(x + h) - y(x), temos a diferença finita progressiva 

de y(x) e a Fórmula Progressiva por: 

y(x) = ~ A y(x) 

h 

Analogamente, podemos obter a Fórmula Regressiva do seguinte modo: 

y'(x)= yix) - yix - h) 
y h 

Considerando Vy(x) = y(x) - y(x - h ), obtemos a diferença finita regressi¬ 
va de y(x) e a Fórmula Regressiva por: 

y(x) - T-Vy(x) 
h 

Podemos também, usando a expansão da função y(x) por Série de Taylor 
nas vizinhanças do ponto x, até ordem 1, desprezando os termos da ordem 
0(h 2 ) e subtraindo y(x - h) de y(x + h), obter: 

, (x) y(x + h)-y(x-h) 
y 2h 

Considerando ôy(x) = y(x + h) - y(x - h) , temos a diferença finita central 
de y(x) e a Fórmula Central dada por: 







/'(*) = 


8y(x) 

2h 


Como anteriormente, podemos obter uma aproximação para y "(x), toman¬ 
do a expansão da função y(x) por Série de Taylor até a ordem 3 e considerando 
a soma das expressões para y(x + h) e y(x - h) , temos: 


y"(x) 


y(x + h)~ 2y(x) + y(x - h) 
h 2 




Observação 

As notações Ay(x), Vy(x) e ôy(x) representam os operadores de diferença 
finita progressivo, regressivo e central, respectivamente, aplicados à função y(x). 

Essas fórmulas de diferenças finitas desenvolvidas para aproximação das 
derivadas de uma função de uma variável podem ser utilizadas para aproxima¬ 
ção das derivadas parciais de uma função de várias variáveis. 

Seja u(x,t) uma função de duas variáveis (x,t), temos as seguintes 
diferenças: 

Fórmula Progressiva na variável t 


^ Ms W**l-»W) = l wo 

K K 


Fórmula Regressiva na variável t 


„„u(x,t)-u(x,t-k)_ 

U t * / — / — 



U(x, t ) 


Fórmula Central na variável x 


U(X + h,t)-U(X-h,t) = ± 
x 2h 2h x 


U,a(X,t) = 


u(x + h,t)- 2u(x,t) + u(x - h,t) 1 2 


h“ 


= -j 1 2 5 x U ( X ’ t) 










6.5.1 Método de Diferenças Finitas 


Seja uma função de duas variáveis u{x,t), considere o plano subdividido 
em retângulos iguais de lados A x = h e A t = k. 

Considere um ponto P = (x, t) na malha, então este será representado por 
P = (x, t) = (x„ tj ) = (ih, jk) . Os índices i e j são inteiros e os valores h e k são os 
espaçamentos da malha nas direções dos eixos x e t. 

Nesse caso, o valor da função u(x,t) no ponto P será denotado por 
u P = u{x j , tj) = u(ih,jk) = li, j. 

O Método das Diferenças Finitas consiste em definir uma malha na região 
de interesse a ser calculada, os valores da equação diferencial parcial e utilizar 
aproximações das derivadas vistas anteriormente nessa equação, num ponto 
genérico P = (x j ,t j ), conforme mostrado na Figura 6.5. 



Figura 6.5 Malha dos pontos no plano. 

O Método de Diferenças Finitas para uma equação do tipo parabólica é 
obtido pelo desenvolvimento da Série de Taylor da função u(x, t ) na variável t 
(tempo) para formar o quociente de diferenças: 


du 

~dt 


(M y ) = 


U(Xj, t J+k ) — U(Xj, ts) kd 2 U , 

- k - + -2W {X ^ 


( 2 ) 


para algum e (t j; t j+1 ) e a Série de Taylor para função u{x,t) na variável 
espacial x: 
























( 3 ) 



em que Ç, e(x,_„ x /+? ). 

Uma vez substituídas as expressões das derivadas (2) e (3) e desprezando 
os termos de ordem O (k) e de 0(h 2 ) na equação diferencial parcial, teremos 
uma equação de diferenças para ser resolvida, escrevendo as derivadas par¬ 


ciais no ponto P = (x.,í y ) como: 



ii _ 011 j -ii 



6.5.2 Método Explícito 


Descreveremos a seguir um método explícito para resolução de uma equa¬ 
ção diferencial parcial, baseado em diferenças finitas. 

As equações diferenciais parciais de segunda ordem classificam-se em 
hiperbólicas, elípticas e parabólicas. Para facilidade de entendimento exibiremos 
um problema unidimensional de difusão dado pela equação diferencial parabóli¬ 
ca para uma função u(x, t ) definida no retângulo [0, L] x [0, T]: 



(4) 


com condições de contorno u(0, t) = u{L, t) = 0 e 0 < t < T, as quais representam 
os valores da função u nas extremidades x = 0 e x = L, as condições iniciais no 
instante t = 0, u(x, 0) = f(x ) , 0<x<L, e a é o coeficiente de difusão. 

O Método das Diferenças Finitas Explícito consiste em determinarmos uma 
aproximação da solução da equação diferencial parcial envolvendo uma função 
u(x, t ) utilizando, para isto, a derivada parcial de u com relação a variável t (tempo) 
pela diferença progressiva e a derivada parcial de u com relação a variável x (espa¬ 
ço) pela diferença central, conforme malha da Figura 6.6: 
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Figura 6.6 Malha do Método Explícito. 


Inicialmente, construímos uma malha no domínio selecionando um número 
inteiro m>0 e definimos um passo na direção da variável x de tamanho h = l/m. 
Em seguida, selecionamos um passo na direção da variável t de tamanho k. 

Os pontos da malha nesse caso são do tipo (x ; , f y ), em que x, = ih para 


i = 0,1,...,m, etj=jk, para j = 0,1,... 


A equação parabólica (4) implica que para todos os pontos interiores à ma¬ 
lha (. Xj,tj ), para todo / = 1,2,..., m-1 e j = 1,2,..., teremos: 


du 

~dt 


(X i ,tj) = a 2 



(X„tj) 


Então, substituindo as expressões (2) e (3) em (4) temos a função u{x,t) 
nos pontos (x^tj), ou seja, u (y = u{x j3 tj) satisfazendo o sistema de equações 
lineares como segue: 


Uj,j + 1 - U iJ =a 2 U J^i- 2U i,j + 


h £ 


(5) 


Para todo i = 1, 2,..., m-1 e j = 1,2,... 

Essa aproximação leva a um erro de truncamento local no ponto (/, j) para 
esta equação de diferenças que é dado por: 




















A equação de difusão discretizada como em (5) pode ser resolvida para 


u ij+1 conhecendo-se os valores u Lj , u j+1j e u j _ 1j , como segue: 



( 6 ) 


Para todo / = 1, 2,..., m -1 e j = 1,2,... 




(7) 


Se conhecermos uma condicão inicial u{x, 0 ) = f{x ), para todo 0 < x < 1, te¬ 
mos que u i0 =f(x i ), para todo i = 0,1,..., m-1 e as condições de contorno nas 
extremidades u(0, t) = 0 e u{1, t) = 0 temos que u 0 j = u m J , respectivamente. 

Dessa forma, podemos determinar todos os elementos da forma u ír ex¬ 
plicitamente conhecendo as aproximações u i0 e continuando o procedimento 
obtemos de maneira similar os valores u j2 ,u i3 ,... 

Exemplo 6.5 

Considere o problema de difusão com valor inicial e condições de contorno 
dado pela equação unidimensional: 



( 8 ) 


dt dx 


Com as condições de contorno u(0,t) = u(1,t) = 0,t>0, e com as condi¬ 
ções iniciais u(x,0) = sen(nx), 0<x<1. 

Resolvendo numericamente e utilizando o método de diferenças progres¬ 
sivas na variável t e centradas em x, considerando os passos na variável x com 
h = 0.2 e com passo na variável t com k = 0.01. Observe que a relação entre os 
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Considerando esses passos temos: 


x j =0.2i, i = 0,1,2,..., 5 
tj = 0.01j, j = 0,1,2,... 


As condições de contorno do problema nos fornecem u 0] = u Xj = 0 e os da¬ 
dos iniciais, u i0 = sen(n(ih )). Representamos os valores das condições iniciais 
na Tabela 6.3: 


Tabela 6.3 


i 

j 

x, = ih = 0.2 i 

U U 

0 

0 

0 

0 

1 

0 

0.2 

0.5877852524 

2 

0 

0.4 

0.9510565165 

3 

0 

0,6 

0.9510565163 

4 

0 

0.8 

0.5877852522 

5 

0 

1 

0 


Aplicando a expressão da fórmula (7): 


u u+ i = 0.25u i+1 j + 0.5u u + 0.25 Uj^j 

ou 

Ut, i+1 = 0-25(u i+1j + 2u i:j + u^j) 

para calcular os valores da solução aproximada nos pontos da (y-t-f)-ésima 
linha no tempo, conhecendo-se apenas os valores da j-ésima linha. 

Assim, obtemos os elementos u jj+1 para j = 1,2,... e os demais elementos 
a partir das condições iniciais e condições de contorno. 

Escolhendo inicialmente j = 0, obtemos os valores aproximados da solu¬ 
ção u 1p u 21 , u 31 e u 41 na primeira linha da malha na variável t, como segue: 

— 0.25(u 2 o + 2u 1 0 + u Q0 ) 

= 0.25(0.9510565165) + 2(0.58778525244 ) + (0) 


= 0.5316567555 














u 21 = 0.25 (u 30 + 2u 20 + u 10 ) = 0.8602387002 


u 31 = 0.25(u 40 + 2u 30 +u 20 ) = 0.8602387005 
u 41 = 0.25{u 5 0 + 2u 4 0 + u 30 ) =0.5316567555 


Da mesma forma, os valores numéricos aproximados de u(x, t) na segun¬ 
da linha da variável t e x, e em todos os pontos do domínio desejado podem ser 
calculados. 


Como a solução exata desse problema é dada pela expressão 

u(x,t) = e n2t sen(nx), podemos comparar os resultados obtidos numericamen¬ 
te com os valores exatos nos pontos (/, j) e obter os erros correspondentes. 


A estabilidade de um método numérico está relacionada ao crescimento 
ou decaimento dos erros decorrentes das operações aritméticas envolvidas na 
resolução do sistema de equações envolvido. 


Dizemos que um método é estável quando esses erros não aumentam 
quando o número de pontos da malha cresce. 

k 1 

Mostra-se que devemos escolher os valores de h e k de modo que — r < —, 

h 2 2 

e já que a = 1 no exemplo dado, escolhemos h = 0,2 e k = 0,01 e obtemos 
1 ^ 

\ = ^r = 0.25 <0.5, de modo que essa condição de estabilidade é satisfeita 
h 

para o método utilizado. 


Caso o valor da relação entre o espaçamento da malha X fosse maior do 
que 0.5, teríamos problemas de estabilidade e buscaríamos outro método para 
resolução. 


6.6 Exercícios 


1. Usando o Método de Euler, resolva as seguintes equações diferenciais 
com valor inicial: 

a) y' = 4-2x, com y(0) = 2. Usando h = 0.01, calcule y(0.05 ). 

b) y' = 1- y/x, com y(2) = 2. Usando h = 0.1, calcule y (2.4). 

c ) y' = -y/x, com y(0 ) = 20. Usando h = 0.2, calcule y(1.0). 




2. Usando o Método de Euler Aperfeiçoado, calcule a solução aproximada 
para as seguintes equações diferenciais com valor inicial: 


a) y' = -x/y, com y(0) = 20. Usando h = 0.2, calcule y (1-2). 

b) y' = y-x, com y(0) = 2. Usando h = 0.1, calcule y(0.9). 


3. Com o Software Numérico - Referência [1] usando o Método de Euler, 
calcule a solução aproximada para as seguintes equações diferenciais com valor 
inicial: 

a) y' = -2y + 1, com y(0) = 1. 

Usando h = 0.1 e depois h = 0.01, calcule y(0.9). O que você pode afir¬ 
mar sobre a qualidade dos resultados obtidos? Justifique suas afirmações 
teoricamente. 

b) y' = y, com y(0) = 1 . 

Usando h = 0.2, calcule y(2.0) e um limitante superior para o erro. 

4. Escreva as equações discretizadas da equação diferencial parabólica da 
difusão para: 

a ) h = 0.1 e k = 0.1 

b ) h =0.5e k= 0.3 

5. Discretize a equação diferencial parcial parabólica da difusão, 


du 

~di 


(x,t) 


d 2 U 

dx 2 


(X,t) 


utilizando o Método das Diferenças Finitas para valores de passos no tempo 
k = 0.05 e no espaço h = 0.1. 

Considere as condições de contorno u(0,t) = u(2,t) = 0, 0 <t <T e as con¬ 
dições iniciais u(x,0) = 1,0, 0<x<2 para calcular o valor de u(x,t) nos pontos 
(0.9, 0.05), (1,0.05), (1.1,0.05) e (1,0.1). 


6. Discretize a equação diferencial parcial parabólica da difusão, 


du 

~di 


(x,t) = a 2 


d 2 u 

dx 2 


(X,t) 
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utilizando o Método das Diferenças Finitas para valores do coeficiente de difusão 
a = 1. 






Escolha valores para os passos no tempo k e no espaço h e encontre um 
valor aproximado para a solução numa região quadriculada de lado 1, conside¬ 
rando as condições de contorno u(0,t) = u(1,t) = 0, 0 <t <T e as condições ini¬ 
ciais u(x,0) = 1,<x<1. Plote os valores das aproximações obtidas nos pontos 
da malha. 

7. Faça um Mapa Conceituai detalhado sobre métodos numéricos para as 
Equações Diferenciais Ordinárias, introduzindo labeis indicando se aprendeu 
(A), não aprendeu (N), se gostou (G), detestou (D ), achou interessante (I), etc. 

8. Faça um Mapa Conceituai detalhado sobre métodos numéricos para as 
Equações Diferenciais Parciais, introduzindo labeis indicando se aprendeu (A), 
não aprendeu (N), se gostou (G), detestou (D ), achou interessante (I), etc. 
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